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ABSTRACT 

We calculate steady-state, one-dimensional hydrodynamic profiles of hot gas in slowly 
accreting (“quiescent”) galactic nuclei for a range of central black hole masses M,, 
parametrized gas heating rates, and observationally-motivated stellar density profiles. 
Mass is supplied to the circumnuclear medium by stellar winds, while energy is injected 
primarily by stellar winds, supernovae, and black hole feedback. Analytic estimates are 
derived for the stagnation radius (where the radial velocity of the gas passes through 
zero) and the large scale gas inflow rate, M, as a function of M. and the gas heating ef¬ 
ficiency, the latter being related to the star-formation history. We assess the conditions 
under which radiative instabilities develop in the hydrostatic region near the stagna¬ 
tion radius, both in the case of a single burst of star formation and for the average 
star formation history predicted by cosmological simulations. By combining a sample 
of measured nuclear X-ray luminosities, Lx, of nearby quiescent galactic nuclei with 
our results for M(M .) we address whether the nuclei are consistent with accreting in 
a steady-state, thermally-stable manner for radiative efficiencies predicted for radia- 
tively inefficiency accretion flows. We find thermally-stable accretion cannot explain 
the short average growth times of low mass black holes in the local Universe, which 
must instead result from gas being fed in from large radii, due either to gas inflows or 
thermal instabilities acting on larger, galactic scales. Our results have implications for 
attempts to constrain the occupation fraction of SMBHs in low mass galaxies using 
the mean Lx — M , correlation, as well as the predicted diversity of the circumnuclear 
densities encountered by relativistic outflows from tidal disruption events. 
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1 INTRODUCTION 

Supermassive black holes (SMBHs) lurk in the centres of 
most, if not all nearby galaxies (see reviews by, e.g. Kor- 
mendy & Richstone 1995; Ferrarese & Ford 2005). However, 
only a few percent of these manifest themselves as luminous 
active galactic nuclei (AGN). Nearly quiescent SMBHs, such 
as those hosting low luminosity AGN, constitute a silent ma¬ 
jority (e.g. Ho 2009). 

Understanding why most SMBHs appear to be inac¬ 
tive requires characterizing their gaseous environments. Gas 
near the SMBH sphere of influence, hereafter denoted the 
‘circumnuclear medium’ (CNM), controls the mass accretion 
rate, M.. The accretion rate in turn determines the SMBH 
luminosity and the feedback of its energy and momentum 
output on larger scales. Dense gas in the nucleus may lead 
to runaway cooling, resulting in bursty episodes of star for¬ 
mation and AGN activity (e.g. Ciotti & Ostriker 2007). 

* E-mail: ag@astro.columbia.edu 


Knowledge of how M. depends on the SMBH mass, M ., 
and other properties of the nucleus informs key questions 
related to the co-evolution of SMBHs and their host galax¬ 
ies with cosmic time (e.g. Kormendy & Ho 2013a; Heck¬ 
man & Best 2014). In the low redshift Universe, SMBH 
growth is dominated by low mass black holes, M. < 10 8 Mq 
(e.g. Heckman et al. 2004), a fact often attributed to the 
trend of ‘cosmic down-sizing’ resulting from hierarchical 
structure growth (e.g, Gallo et al. 2008). However, the phys¬ 
ical processes by which typical low mass black holes ac¬ 
crete could in principle be distinct from those operating at 
higher SMBH masses, or those in AGN. Of key importance 
is whether SMBHs grow primarily by the accretion of gas fed 
in directly from galactic or extragalactic scales, or whether 
significant growth can result also from local stellar mass loss 
in the nuclear region. 

A better understanding of what mechanisms regulate 
accretion onto quiescent SMBHs would shed new light on a 
variety of observations, such as the occupation fraction of 
SMBHs in low mass galaxies. Miller et al. (2015) use the av- 


© 0000 The Authors 


2 Generozov, Stone, & Metzger 


erage relationship between the nuclear X-ray luminosities, 
Lx, of a sample of early type galaxies and their associated 
SMBH masses to tentatively infer that the SMBH occupa¬ 
tion fraction becomes less than unity for galaxies with stellar 
masses M* < 10 1o Mq (M m < 10' Mq). This method relies 
on extrapolating a power-law fit of the Lx — M. distribu¬ 
tion to low values of Lx below the instrument detection 
threshold, an assumption which is questionable if different 
physical processes control the accretion rates onto the lowest 
mass SMBHs. 

The gas density in galactic nuclei also influences the 
emission from stellar tidal disruption events (TDEs), such 
as the high energy transient Swift J1644+57 (Bloom et al. 
2011, Burrows et al. 2011; Levan et al. 2011; Zauderer et al. 
2011). This event was powered by an impulsive relativistic 
jet, which produced synchrotron radio emission as the jet 
material decelerated from shock interaction with the CNM 
of the previously quiescent SMBH (Giannios & Metzger 
2011; Zauderer et al. 2011). Modeling of J1644+57 showed 
that the CNM density was much lower than that measured 
surrounding Sgr A* on a similar radial scale (Berger et al. 
2012; Metzger et al. 2012). However, a TDE jet which en¬ 
counters a denser CNM would be decelerated more rapidly, 
producing different radio emission than in J1644+57. Vari¬ 
ations in the properties of the CNM could help explain why 
most TDEs appear to be radio quiet (e.g. Bower et al. 2013; 
van Velzen et al. 2013). 

Gas comprising the CNM of quiescent (non-AGN) 
galaxies can in principle originate from several sources: (1) 
wind mass loss from predominantly evolved stars; (2) stellar 
binary collisions; and (3) unbound debris from a recent TDE. 
Stellar wind mass loss is probably the dominant source inso¬ 
far as collisions are relevant only in extremely dense stellar 
environments for very young stellar populations (Rubin & 
Loeb 2011), while MacLeod et al. (2013) find that TDEs con¬ 
tribute subdominantly to the time-averaged accretion rate 
of quiescent SMBHs. 

The gas inflow rate on large scales is much easier to con¬ 
strain both observationally and theoretically than the black 
hole (horizon scale) accretion rate. Ho (2009) determines the 
inflow rates in a sample of early-type galaxies by using X- 
ray observations to determine the Bondi accretion rate, and 
also by using estimated mass loss rates of evolved stars. Both 
methods lead him to conclude that the available gas reservoir 
is more than sufficient to power the observed low-luminosity 
AGN, assuming the standard ~ 10 per cent radiative effi¬ 
ciency for thin disc accretion. Several lines of evidence now 
suggest that low-luminosity AGN result from accretion pro¬ 
ceeding in a radiatively inefficient mode (Yuan & Narayan 
2014), due either to the advection of gravitationally-released 
energy across the SMBH horizon (e.g. Narayan & Yi 1995) or 
due to disc outflows, which reduce the efficiency with which 
the inflowing gas ultimately reaches the SMBH (e.g. Bland- 
ford & Begelman 1999; Li et al. 2013). 

Another approach to determine the inflow rates, which 
we adopt in this paper, is to directly calculate the density, 
velocity and temperature profiles of the CNM using a phys¬ 
ically motivated hydrodynamic model. Mass is injected into 
the nuclear environment via stellar winds, while energy is in¬ 
put from several sources including stellar winds, supernovae 
(SNe), and AGN feedback (Quataert 2004; De Colle et al. 
2012; Shcherbakov et al. 2014). Unlike previous works, which 


focused primarily on modeling individual galaxies, here we 
model the CNM properties across a representative range of 
galaxy properties, including different SMBH masses, stellar 
density profiles, and star formation histories (SFHs). 

Previous studies, employing multi-dimensional numeri¬ 
cal hydrodynamics and including variety of (parametrized) 
physical effects, have focused on massive elliptical galax¬ 
ies (e.g. Ciotti & Ostriker 2007; Ciotti et al. 2010). These 
works show the periodic development of cooling instabilities 
on galactic scales, which temporarily increase the gas in¬ 
flow rate towards the nucleus until feedback becomes strong 
enough to shut off the flow and halt SMBH growth. 

In this paper we focus on time-independent models, in 
which the nuclear gas receives sufficient heating to render 
radiative cooling negligible. This approach allows us to sys¬ 
tematically explore the relevant parameter space and to de¬ 
rive analytic expressions that prove useful in determining 
under what conditions cooling instabilities manifest in the 
nuclear region across the expected range of galaxy prop¬ 
erties, or whether other (non-AGN) forms of feedback can 
produce a prolonged state of steady, thermally stable accre¬ 
tion. Even if cooling instabilities develop on galactic scales 
over longer ~ Gyr time-scales, we aim to explore whether a 
quasi steady-state may exist between these inflow events on 
smaller radial scales comparable to the sphere of influence. 

In the presence of strong heating, one-dimensional 
steady-state flow is characterized by an inflow-outflow struc¬ 
ture, with a critical radius known as the “stagnation radius” 
r s where the radial velocity passes through zero. Mass loss 
from stars interior to the stagnation radius is accreted, while 
that outside r a is unbound in an outflow from the nucleus. 
The stagnation radius, rather than the Bondi radius, thus 
controls the inflow rate (although we will show that r a usu¬ 
ally resides near the nominal Bondi radius). When heating is 
sufficiently weak, however, the stagnation radius may move 
to much larger radii or not exist at all, significantly increas¬ 
ing the inflow rate the SMBH, i.e. a “cooling flow”. However, 
the hydrostatic nature of gas near the stagnation radius also 
renders the CNM at this location particularly susceptible to 
local thermal instabilities, the outcome of which could well 
be distinct from the development of a cooling flow. 

This paper is organized as follows. In §2 we describe 
our model, including the sample of galaxy properties used 
in our analysis (§2.1) and our numerical procedure for cal¬ 
culating the steady-state hydrodynamic profile of the CNM 
(§2.2). In §3 we describe our analytic results, which are jus¬ 
tified via the numerical solutions we present in §4. We move 
from a general and parametrized treatment of heating to a 
physically motivated one in §5, where we consider a range 
of physical processes that can inject energy into the CNM. 
In §6 we discuss the implications of our results for topics 
which include the nuclear X-ray luminosities of quiescent 
black holes, jetted TDEs, and the growth of SMBHs in the 
local Universe. In §7 we summarize our conclusions. Table 
1 provides the definitions of commonly used variables. Ap¬ 
pendix A provides useful analytic results for the stagnation 
radius, while Appendix B provides the details of our method 
for calculating stellar wind heating and mass input. 
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Table 1. Definitions of commonly used variables 


Variable Definition 


M. 

Vw 

V%u 

^0 

^dyn 

% 

C 

r s 

^inf 

r b 

Ha 

P*{r) 

p(r) 

M+(r) 

M eac (y ) 
q(r) 

V 

T* 

r 

s 

M 

M. 

Mi. 

Me 

JVf'Pi 

9heat /1 Qrad I 


Black hole mass 

Total heating parameter, including minimum heating rate from stellar and black hole velocity dispersion 

Total heating parameter, excluding minimum heating rate from velocity dispersion 

Stellar velocity dispersion (assumed to be radially constant). 

t / ctq : Dynamical time at large radii (where stellar potential dominates) 

(r 3 /(GM#)) 1 / 2 : Free fall time (where the black hole potential dominates) 

Alternative heating parameter, £ = y/l + ( v w /(jq ) 2 (eq. [16]) 

Stagnation radius, where gas radial velocity goes to zero 
Radius of sphere of influence (eq. [3]) 

Outer break radius of stellar density profile 

Radius interior to which SN la are infrequent compared to the dynamical time-scale (eq. [36]) 

3D radial stellar density profile 
Gas density of CNM 

Total enclosed stellar mass inside radius r 

Total enclosed mass inside radius r (SMBH + stars) 

Mass source term due to stellar winds, q oc p* (eq. [8]) 

Parameter setting normalization of mass input from stellar winds (eq. [8]) 

Age of stellar population, in case of a single burst of star formation 

Power-law slope of radial stellar surface brightness profile interior to the break radius 

Power-law slope of the 3D stellar density profile inside of the break radius, <5 = T + 1. 

Large scale inflow rate (not necessarily equal to the SMBH accretion rate) 

SMBH accretion rate, M. = qM, where a < 1 accounts for outflows from the accretion disc on small scales. 
Maximum accretion rate as limited by SN la (eq. [39]) 

Equilibrium inflow rate set by Compton heating acting alone (eq. [46]) 

Maximum accretion rate for thermally stable accretion (eq. [29]) 

Ratio of external heating (stellar winds, SN la, MSPs) to radiative cooling (eq. [25]) 

Hubble time-scale 

Gas density power-law slope at the stagnation radius (eq. [13]) 


2 MODEL 

2.1 Galaxy models 

Lauer et al. (2007) use Hubble Space Telescope WFPC2 
imaging to measure the radial surface brightness profiles for 
hundreds of nearby early type galaxies. The measured profile 
is well fit by a “Nuker” law parameterization: 

m = h2v-™°c T a +n- (/J - r)/a , * = « 

' 1, 

i.e., a broken power law that transitions from an inner power 
law slope, r, to an outer power law slope, /3, at a break 
radius, ?- b . If the stellar population is spherically symmetric, 
then this corresponds to a 3D stellar density p* oc r~ 1_r for 
r Tb and p* oc r -1 ~' 3 for r rb- We can write the 
deprojected stellar density approximately (formally, this is 
the a —> oo limit) as 

„ _ fp*ln=f (r/rinf)“ 1_r r<r b 
P* — ) | / / \-i-p . (2) 

yr/r h ) p r > r b , 

where p*|r inf is the stellar density at the radius of the black 
hole sphere of influence 1 , 

nnf =i GM./a 2 . « 14M. 0 ;! pc, (3) 


where M.,8 = M,/10 s Mq and the second equality in (3) 
employs the M, — a relationship of McConnell et al. (2011), 


M. ~ 2 x 10 8 


200 km s - 




(4) 


This may be of questionable validity for low mass black holes 
(e.g., Greene et al. 2010; Kormendy & Ho 2013b). Also, 
several of the black hole masses used in McConnell et al. 
(2011) were underestimated (Kormendy & Ho 2013b). How¬ 
ever, our results are not overly sensitive to the exact form 
of the M, — a relationship that we use. 

A galaxy model is fully specified by four parameters: 
M m , r, i'b, and /3. We compute models for three different 
black hole masses, Al, = 10 6 , 10 7 , 10 8 M@. The distribution 
of F in the Lauer et al. (2007) sample is bimodal, with a 
concentration of “core” galaxies with T < 0.3 and a concen¬ 
tration of “cusp” galaxies with V > 0.5. We bracket these 
possibilities by considering models with F = 0.8 and T = 0.1. 

We fix j3 = 2 but find that our results are not overly 
sensitive to the properties of the gas flow on radial scales > 
rb- The presence of the break radius is, however, necessary 
to obtain a converged steady state for some regions of our 
parameter space 2 . We consider solutions calculated for up 
to four values of r b : 50 pc, 100 pc, 200 pc, and 400 pc, 


1 This approximate expression agrees surprisingly well with the 
mean empirical scaling relation for r; n f(M.) calculated in ap¬ 
pendix C of Stone & Metzger (2014), although we note that this 
relation has significant scatter. Also, we note that the scaling 
is somewhat different for cores and cusps. Fixing the power law 


slope, r'; n f ~ 25(8 )Mfg for core (cusp) galaxies. We will use sepa¬ 
rate scaling relations for r , n f for cores and cusps unless otherwise 
noted. 

2 In addition, for /3 < 2 the stellar density must steepen at still 
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motivated by the range of break radii from the Lauer et al. 
(2007) sample . 3 We neglect values of rb which would give 
unphysically large bulges for a given M,. 

Finally, we note that Lauer et al. 2005 find that ~ 60% 
of cusp galaxies and ~ 29% of core galaxies have (gener¬ 
ally unresolved) emission in excess of the inward extrap¬ 
olation of the Nuker law. Indeed, some low mass galaxies 
with M. < 10 8 Mg possess nuclear star clusters (Graham 
& Spitler 2009), which are not accounted for by our simple 
parametrization of the stellar density. In such cases gas and 
energy injection could be dominated by the cluster itself, i.e. 
concentrated within its own pc-scale “break radius” which is 
much smaller than the outer break in the older stellar pop¬ 
ulation on much larger scales. Although our analysis does 
not account for such an inner break, we note that for high 
heating rates the stagnation radius and concomitant inflow 
rate are not sensitive to the break radius. 


2.2 Hydrodynamic Equations 

Following Quataert (2004, see also Holzer & Axford 1970; 
De Colie et al. 2012; Shcherbakov et al. 2014), we calculate 
the density p, temperature T, and radial velocity v of the 
CNM for each galaxy model by solving the equations of one¬ 
dimensional, time-dependent hydrodynamics, 


dp Id, 2 \ 
dt+^dr^ V ^ =q 
dv dv \ dp 
dt. +V dr) ~ 


'ds ds\ 
pT [dt +V dr) =q 


GM el 


: P~ 

rp2 

- qv 

v 2 

J + 

C4 3 

7 V 

2 

7-1 P. 


(5) 

( 6 ) 
(7) 


where p and s are the pressure and specific entropy, 
respectively, and M enc = M+(r) + M. is the enclosed mass 
(we neglect dark matter contributions). We adopt an ideal 
gas equation of state with p = pkT / pm p with /t = 0.62 and 
7 = 5/3. The source term in equation (5), 


q = 


VP* 
th ’ 


( 8 ) 


represents mass input from stellar winds, which we 
parametrize in terms of the fraction p of the stellar density 
p» being recycled into gas on the Hubble time th = 1.4 x 10 10 
yr. To good approximation p ~ 0.02(r*/th ) -1,3 at time t* 
following an impulsive starburst (e.g., Ciotti et al. 1991) 4 , 
although p is significantly higher for continuous SFHs (bot¬ 
tom panel of Fig. 7). 

Source terms oc q also appear in the momentum and 
entropy equations (eqs. [ 6 ] and [7]) because the isotropic in¬ 
jection of mass represents, in the SMBH rest frame, a source 


larger radii to avoid the mass enclosed diverging to infinity. How¬ 
ever, we do not find it necessary to include this outer break. 

3 For the core galaxies the break radius follows the scaling rela¬ 
tionship rb ~ 90 (A/, g) 0 ' 5 pc, with scatter of approximately one 
dex. Most of the cusp galaxies have rb between 100 and 1000 pc, 
and lack a clear trend with M ,. The mean rb for cusp galaxies in 
the Lauer et al. (2007) sample is ~ 240 pc. 

4 Ciotti et al. (1991) give the mass return rate from evolved stars 
as a function of B-band luminosity instead of volumetrically, but 
our expressions are equivalent. 


of momentum and energy relative to the mean flow. Physi¬ 
cally, these result from the mismatch between the properties 
of virialized gas injected by stellar winds and the mean back¬ 
ground flow. The term tx p/p = c 2 is important because it 
acts to stabilize the flow against runaway cooling (§3.4). 

The term oc v%, = <r(r) 2 -|-i>2, in the entropy equation ac¬ 
counts for external heating sources (e.g., Shcherbakov et al. 
2014), where 


3 GM. 9 

+ °o> 


(F + 2 )r 


(9) 


is the stellar velocity dispersion. This accounts for the min¬ 
imal amount of shock heating from stellar winds due to the 
random motion of stars in the SMBH potential. We take cro 
to be constant and use a 2 « 3er^, where a. ~ 170M^'f km/s 
is the velocity dispersion from the McConnell et al. (2011) 
M, — a relation. The second term, v%,, parametrizes addi¬ 
tional sources of energy input, including faster winds from 
young stars, millisecond pulsars (MSPs), supernovae, AGN 
feedback, etc (§5). We assume that v w is constant with ra¬ 
dius, i.e. that the volumetric heating rate is proportional to 
the local stellar density. 

Our model does not take into account complications 
such as more complicated geometries or the discrete na¬ 
ture of real stars (Cuadra et al. 2006, 2008). In the case 
of Sgr A* these effects reduce the time-averaged inflow rate 
to ~ lO _ 6 M 0 yr _1 - an order of magnitude less than a ID 
spherical model (Cuadra et al. 2006). Motions of individual 
stars can produce order of magnitude spikes in the accretion 
rate (Cuadra et al. 2008). 

To isolate the physics of interest, our baseline calcula¬ 
tions neglects three potentially important effects: heat con¬ 
duction, radiative cooling, and rotation. Heat conduction 
results in an an additional heating term in equation (7), 


9cond = V • (kVT), (10) 

where k = K sp i tz /(l + VO (Dalton & Balbus 1993) is the 
conductivity and K; sp itz = kqT 5 ' 2 is the classical Spitzer 
(1962) value (ko — 2 x 10 -6 in cgs units). The flux limiter 
V> = K S pitzVT/(5</pc 3 ) saturates the conductive flux if the 
mean free path for electron coulomb scattering exceeds the 
temperature length scale, where c s = (kT / pm v ') 1 ^ 2 is the 
isothermal sound speed and </ < 1 is an uncertain dimen¬ 
sionless constant (we adopt 0 = 0.1). Even a weak magnetic 
field that is oriented perpendicular to the flow could sup¬ 
press the conductivity by reducing the electron mean free 
path. However, for radially-decreasing temperature profiles 
of interest, the flow is susceptible to the magneto-thermal 
instability (Balbus 2001), the non-linear evolution of which 
results in a radially-directed field geometry (Parrish & Stone 
2007). In §3.3 we show that neglecting conductivity results 
in at most order-unity errors in the key properties of the 
solutions. 

Radiative cooling contributes an additional term to 
equation (7), of the form 

qr ad = —A(T)n 2 , (11) 

where n = p/pm p and A(T) is the cooling function. We ne¬ 
glect radiative cooling in our baseline calculations, despite 
the fact that this is not justified when the wind heating 
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v w is low or if the mass return rate 77 is high. Once radia¬ 
tive cooling becomes comparable to other sources of heat¬ 
ing and cooling, its presence can lead to thermal instabil¬ 
ity (e.g. Gaspari et al. 2012b, McCourt et al. 2012, Li & 
Bryan 2014a) that cannot be accounted for by our ID time- 
independent model. Our goal is to use solutions which ne¬ 
glect radiation to determine over what range of conditions 
cooling instabilities will develop (§3.4). 

Equations (5)-(7) are solved using a sixth order finite 
difference scheme with a third order Runge-Kutta scheme for 
time integration and artificial viscosity terms in the velocity 
and entropy equations for numerical stability (Brandenburg 
2003) 5 . We assume different choices of v w = 300,600,1200 
km s -1 spanning a physically plausible range of thermally 
stable heating rates. Although we are interested in the 
steady-state inflow/outflow solution (assuming one exists), 
we solve the time-dependent equations to avoid numerical 
issues that arise near the critical sonic points. 

Our solutions can be scaled to any value of the mass 
input parameter, 77 , since the mass and energy source terms 
scale linearly with p or p*; however, the precise value of 77 
must be specified when cooling or thermal conduction are 
included. We check the accuracy of our numerical solutions 
by confirming that mass is conserved across the grid, in ad¬ 
dition to the integral constraint on the energy (Bernoulli 
integral). 


3 ANALYTIC RESULTS 

We first describe analytic estimates of physical quantities, 
such as the stagnation radius r s and the mass inflow rate, 
the detailed derivation of which are given in Appendix A. 


3.1 Flow Properties Near the Stagnation Radius 


Continuity of the entropy derivative at the stagnation radius 
where v = 0 requires that the temperature at this location 
be given by (eq. [A2]) 


T\r„ 


7 — 1 prripV 2 ^ 7 — 1 p,m v Vw 13 + 8 r 

7 2k v m ^>a 0 7 2k 13 + 8 L — 6 v 


5.0 x 10 6 U 500 K core 
5.5 x 10 6 U 500 K cusp, 


( 12 ) 


where U500 = v w /(500 km s -1 ) and v = —d\np/d\nx\ re is the 
density power-law slope at r = r s . Empirically, we find from 
our numerical solutions that 

^~i(4r + 3) (13) 

Hydrostatic equilibrium likewise determines the value 
of the stagnation radius (Appendix A, eq. A 6 ) 


r s 


GM, [~ 13 + sr 
uvl, 4 + 2r 


GM. 

v(vl + erg) 


' M*L S 13 + 8r 

4-— -I- 

M. 4 + 2r 


3 v 

2 + r 


(14) 


5 code is available at https://github.com/alekseygenerozov/hydro 


For high heating rates v w 00 , the stagnation radius re¬ 
sides well inside the SMBH sphere of influence. In this case 
M i ,\ rs /M. <C 1, such that equation (14) simplifies to 

/ 13 + 8 r 3 v \ GM. 

’'raijo'o V 4 + 2r 2 + ry wZ, 

f8pc M.^v^yr" 1 core 
|4pc M., 8 v^ 0 2 0 yr _1 cusp, 

where we have used equation (13) to estimate v separately 
for core (r = 0 . 1 ; v « 1 ) and cusp (r = 0 . 8 ; v « 0 . 6 ) galax¬ 
ies. This expression is similar to that obtained by Volonteri 
et al. ( 2011 ) on more heuristic grounds (their eq. 6 ). 

In the opposite limit of weak heating (v w < no) the stag¬ 
nation radius moves to large radii, approaching the break 
radius rb in the stellar density profile, implying that all of 
the interstellar medium (ISM) inside of rb is inflowing. In 
particular, we find that r s approaches rb for heating below 
a critical threshold 



(16) 


Condition (16) approximately corresponds to the re¬ 
quirement that the heating rate exceed the local escape 
speed at the break radius, rb. This result makes intuitive 
sense: gas is supplied to the nucleus by stars which are grav¬ 
itationally bound to the black hole, so outflows are possible 
only if the specific heating rate ~ ug, significantly exceeds 
the specific gravitational binding energy. 


3.2 Inflow Rate 

The large scale inflow rate towards the SMBH is given by 
the total mass loss rate interior to the stagnation radius 
(eq. [ 8 ]), 

M = 4 - 7 T P qr 2 dr = = ^.(r s /r inf ) 2 - r 

Jo 

4.5 x 10 _ 5 M. 1 ;J 6 u^ 0 3 0 ' 8 ? 7 o.o 2 M 0 yr " 1 core 
3.2 x 10 _ 5 Mg;| 8 u^ 0 o 4 ? 7 o.o 2 M 0 yr _1 cusp, 

where we have assumed v m 00 by adopting equation (15) 
for r s . The resulting Eddington ratio is given by 

M _ f 2.0 x 10 - 5 M. 0 ; 8 76 u 5 - 0 3 0 8 77 o.o 2 

Me dd ~ \l.4x 10- 5 M.°; g 48 r 5 - 0 2 0 Vo2 

where M e dd = 2.2 M. j 8 M 0 yr _1 is the Eddington accretion 
rate, assuming a radiative efficiency of ten per cent. Note the 
sensitive dependence of the inflow rate on the wind heating 
rate. Equation (18) is the radial mass inflow rate on rela¬ 
tively large scales and does not account for outflows from 
the SMBH accretion disc (e.g. Blandford & Begelman 1999; 
Li et al. 2013), which may significantly reduce the fraction 
of M that actually reaches the SMBH. Thus we distinguish 
between the large scale inflow rate, M, and the accretion 
rate onto the black hole, M.. 

The gas density at the stagnation radius, p\ r ., is more 
challenging to estimate accurately. By using an alternative 
estimate of M as the gaseous mass within the stagnation 


cusp, 


(18) 
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radius divided by the free-fall time tg| rB = (Vg/GM .) 1 ^ 2 
(47r/3)r 3 p| rs 


M 


tff|r 


in conjunction with eqs. (17) and (15), we find that 


P\r 


5.2 x 10 26 M. g' 2 u 50 g 8 po.o 2 gem 3 core, 


1.0 x 10 25 M.g 5 v 500 V 0.02 gem 3 


cusp 


(19) 


( 20 ) 


It is useful to compare our expression for M (eq. [17]) 
to the standard Bondi rate for accretion onto a point source 
from an external medium of specified density and tempera¬ 
ture (Bondi 1952): 

M b = 47rAr|p| T . B Uff| r . B , (21) 


where vb = GAf/c 2 ad is the Bondi radius, c s , a d = 
\J •ykT / pm p is the adiabatic sound speed, Uff| rB = 
vb /Iff|r B = (GA/./rn) 1 ^ 2 and A is a parameter of order 
unity. 

Equation (19) closely resembles the Bondi formula 
(eq. [21]) provided that vb is replaced by r s . Indeed, for 
r s < Tint we have that (eq. [14]) 


13 + sr GM. 13 + sr 
4 + 2r vvi (2 + r)(3 + 4r) 


( 22 ) 


where the second equality makes use of equation (12). 


3.3 Heat Conduction 


Our analytic derivations neglect the effects of heat conduc¬ 
tion, an assumption we now check. The ratio of the magni¬ 
tude of the conductive heating rate (eq. [10]) to the external 
heating rate at the stagnation radius is given by 


V ■ (kVT) 2t h K. 0 T\l / e 2 / K 0 T 7/2 \ r \ 

qvl/2 rs r%rip*\r s vl V 5 r s cppc 3 s J 


{ 20 ^0. l 2 M, rvlol unsaturated (core) 
30T]QQ 2 M~g ,5 U 5 oo unsaturated (cusp) 
2 (j) saturated, 

(23) 

where the second equality makes use of equations (12), (15), 
and we have made the approximations V 2 ~ l/r B , V ~ 1 /r s . 
The stellar density profile is approximated as (eq. [2]) 

m .(2 — r) 

47rr ? n f \n*fj 

f 7.9 x 10 _19 M.J- 2 u|6 2 ogcm -3 core 
1^2.3 x 10 _18 A/“g' 5 u 3 oo gem -3 cusp/ 



where the stagnation radius is assume to reside well inside 
the Nuker break radius. 

Equation (23) shows that, even when conduction is sat¬ 
urated, our neglecting of heat conduction near the stagna¬ 
tion radius results in at most an order unity correction for 
causal values of the saturation parameter 0 < 0.1. Our nu¬ 
merical experiments which include conductive heating con¬ 
firm this (§4). We do not consider the possibility that the 
conduction of heat from the inner accretion flow can affect 
the flow on much larger scales (Johnson & Quataert 2007). 



Figure 1. Minimum effective wind heating parameter required 
for thermal stability as a function of SMBH mass. Black lines show 
uti (eq. [27]), the heating rate required for (9heat/l<jradl)r a > 10 
in the high-heating limit when the stagnation radius lies interior 
to the influence radius, for different values of the mass loss param¬ 
eter r; as marked. Blue lines show the minimum heating parameter 
required to have C > Cc = (ri,/r; n f) 0 - 5 ( 1—r l (eq. [16]), separately 
for cusp (solid) and core (dashed) galaxies. Based on the Lauer 
et al. (2007) sample we take rj n f = 25(8 )M® e pc and ry = 90M^ | 
pc (240 pc) for cores (cusps). For ( < £ c the stagnation radius 
moves from inside the influence radius, out to the stellar break 
radius ry,. This renders the flow susceptible to thermal runaway, 
even if v w > uxi- 


3.4 Thermal Instability 

Radiative cooling usually has its greatest impact near or ex¬ 
ternal to the stagnation radius, where the gas resides in near 
hydrostatic balance. If radiative cooling becomes important, 
it can qualitatively alter key features of the accretion flow. 
Initially hydrostatic gas is thermally unstable if the cooling 
time is much less than the free-fall time , potentially result¬ 
ing in the formation of a multi-phase medium (Gaspari et al. 
2012b, 2013, 2015; Li & Bryan 2014b). 

Even if the hot plasma of the CNM does not condense 
into cold clouds, the loss of pressure can temporarily increase 
the inflow and accretion rates by producing a large-scale 
cooling flow. When coupled to feed-back processes which 
result from such enhanced accretion, this can lead to time- 
dependent limit cycle behavior (e.g. Ciotti & Ostriker 2007; 
Ciotti et al. 2010; Yuan & Li 2011, Gan et al. 2014), which 
is also inconsistent with our assumption of a steady, single¬ 
phase flow. 

Cooling instability can, however, be prevented if desta¬ 
bilizing radiative cooling (q ra d oc T~ 2 ' 7 ) is overwhelmed 
by other sources of cooling, namely the stabilizing term 
oc — qc 2 oc — T in the entropy equation (eq. [7]). Neglect¬ 
ing radiative cooling to first order, this term is balanced 
at the stagnation radius by the external heating term, 
(jheat = qVw/ 2- One can therefore assess thermal stability 
by comparing the ratio of external heating to radiative cool- 
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ing |<?rad| (eq. [ 11 ]), 

gheat ^ f 630r)ol 2 M~g 76v lo O , core 

19rad| rs [660?;,);J 2 M“g- 48 tl56o , CUSp, 


(25) 


where we have used equations (20) and (24) for the gas and 
stellar densities at the stagnation radius, respectively. We 
have approximated the cooling function for T < 2 x 10' 
K as A (T) = 1.1 x 1CT 22 (T/10 e K) “°' 7 erg cm 3 s~\ which 
assumes solar metallicity gas (Draine 2011; his Fig. 34.1). 

To within a constant of order unity, equation (25) 
also equals the ratio of the gas cooling time-scale t coo i = 
(3nfcT/2^t)/|q ra d| to the free-fall time iff at the stagnation 
radius. This equivalence can be derived using the equality 


I _ 3gtflf |i~s 

Pks — 2 — r ' 


(26) 


that results by combining equations (17),(19), and (24). 
Sharma et al. (2012) argue cooling instability develops in 
a initially hydrostatic atmosphere if i coo i < lOfft, so equa¬ 
tion (25) represents a good proxy for instability in this case 
as well. 

Based on our numerical results (§4) and the work of 
Sharma et al. (2012) we define thermally stable flows ac¬ 
cording to the criterion q heat > 10 |g rad | (t coD i > 10 t ff ) being 
satisfied near the stagnation radius. This condition trans¬ 
lates into a critical minimum heating rate 


280?7o;o|A7°'g 1 kms 1 ,core 
240??S;o|A/°j( 8 kms -1 ,cusp 


(27) 


Equations (25) and (27) are derived using expressions for 
the stagnation radius and gas density in the high heating 
limit of C > Cc = (Wrw) 0 - 5 ^ 1 ’) (eq. [16]). However, for 
£ < £ c , the stagnation radius diverges to the break radius 
TV The quasi-hydrostatic structure that results in this case 
greatly increases the gas density, which in practice renders 
the flow susceptible to thermal runaway, even if v w > Vti ac¬ 
cording to equation (27). In other words, the true condition 
for thermal instability can be written 


IV, > maxfwxf, oo^b/rinf)*- 1 n ] (28) 

Figure 1 shows vti(M.) for different mass input pa¬ 
rameters 77 , as well as the [ 77 -independent] £ > criterion, 
shown separately for cusp and core galaxies. Based on the 
Lauer et al. (2007) sample we take n n f = 25(8)A/ o e pc and 
r b = 90 Mff (240 pc) and thus r b /ri n f ~ 4M~g 1 (30M~g' e ) 
for cores (cusps). We see that the for high 77 and low M. 
the v w > vti criterion is more stringent, while for low 77 and 
high M,, the C > criterion is more stringent. 

The minimum heating rate for thermal stability corre¬ 
sponds to the maximum thermally-stable inflow rate. From 
equations (27) and (28) this is 


Mti 

ATedd 


— 


2 x 10- 4 77o.olA< 8 36 
mm < _ ’ „ 

2 x 10 °vo.o 2 M .; 8 

8 x 10- 5 t7 
mm < r ’ , 

2 X 10 7?0.02M«'g 


(29) 


cusps. 


Note that since the SMBH accretion rate cannot exceed 
the large scale inflow rate, ALti also represents the maximum 
thermally stable accretion rate. 


What we describe above as “thermal instability” may in 
practice simply indicate an abrupt transition from a steady 
inflow-outflow solution to a global cooling flow, as opposed 
a true thermal instability. In the former case the stagnation 
radius diverges to large radii, increasing the density in the 
inner parts of the flow, which increases cooling and creates a 
large inflow of cold gas towards the nucleus. A true thermal 
instability would likely result in a portion of the hot ISM 
condensing into cold clouds, a situation which may or may 
not be present in a cooling flow. In this paper we do not 
distinguish between these possibilities, although both are 
likely present at some level. Finally, note that if the CNM 
were to “regulate” itself to a state of local marginal thermal 
instability (as has previously been invoked on cluster scales; 
e.g., Voit et al. 2015), then equation (29) might naturally 
reflect the characteristic mass fall-out rate and concomitant 
star formation rate. 


3.5 Angular Momentum 

Our spherically symmetric model neglects the effects of an¬ 
gular momentum on the gas evolution. However, all galaxies 
possess some net rotation, resulting in centrifugal forces be¬ 
coming important at some radius r c i rc = ll/(GM,). Here 
l s = ( rV < j,)\r s is the stellar specific angular momentum near 
the stagnation radius, from which most of the accreted mass 
originates, where V), is the stellar azimuthal velocity. 

Emsellem et al. (2007) use two-dimensional kinematic 
data to measure the ratio of ordered to random motion in a 
sample of early type galaxies, which they quantify at each 
galactic radius R by the parameter 

A 8 ,A ~ \ (30) 

{RVV 2 + a 2 ) cr V ’ 

where a is the velocity dispersion and the brackets indicate 
a luminosity-weighted average. The circularization radius of 
the accretion flow can be written in terms of Xr as 


T'circ r s 2 / ,2 /ni \ 

-lA 1 ! 

7* s 7-inf 

where we have used the definition r bl f = GM,/al and in 
the second inequality have assumed that r s < ri n f, a condi¬ 
tion which is satisfied for the thermally-stable solutions of 
interest. 

Emsellem et al. (2007) (their Fig. 2) find that A r is 
generally < 0.1 on radial scales < 10 per cent of the galaxy 
half-light radius and that Xr decreases with decreasing R 
interior to this point. From equation (31) we thus conclude 
that r c Re < O.Olrv. For the low inflow rates considered the 
gas ( M/M^dd < 0 . 01 ) would be unable to cool on a dy¬ 
namical time and would likely drive equatorial and polar 
outflows (Li et al. 2013). Our model cannot capture such 
two dimensional structures, but would still be relevant for 
intermediate polar angles where the gas is inflowing. 


4 NUMERICAL RESULTS 

Our numerical results, summarized in Table 2, allow us to 
study a range of CNM properties and to assess the validity 
of the analytic estimates from the previous section. 

Figure 2 shows profiles of the density p(r), temperature 
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Table 2. Summary of Numerical Solutions 


M. 

(M 0 ) 


Vw 

(km s _1 ) 

r (“) 

'b 

(pc) 

r /J b '> 

rs / r inf 



Unstable 77’s 

U edd i r , ■ ’ 

V 14rad 1 / Vs [V U ^ j 

Cusp Galaxies, F = 

0.8 







10 6 


300 

- 

0.12 

5.1 x 10 -5 

28 

0.6 



600 

- 

3.2 x 10“ 2 

1.0 x 10“ 5 

1.9 x 10 3 




1200 

- 

7.9 x 10 -3 

1.8 x 10 -6 

7.7 x 10 4 


10 7 


300 

50 

0.38 

2.0 x 10“ 4 

4.5 

0.2,0.6 




100 

TI 

TI 

TI 





200 

TI 

TI 

TI 





400 

TI 

TI 

TI 




600 

25 

8.1 x 10“ 2 

3.2 x 10 -5 

6.2 x 10 2 





100 

8.1 x 10 -2 

3.2 x 10 -5 

6.2 x 10 2 




1200 

100 

2.0 x 10“ 2 

6.0 x 10 -6 

2.6 x 10 4 


10 8 


300 

25 

0.69 

4.2 x 10 -4 

5.5 

0.6 




50 

0.92 

6.0 x 10“ 4 

2 

0.6 




100 

1.4 

9.5 x 10 -4 

0.48 

0.2,0.6 




200 

2.5 

1.9 x 10 -3 

5.3 x 10 -2 

0.2,0.6 



450 

100 

0.43 

2.4 x 10 -4 

19 


...t 


450 

100 

1.2 

8.1 x 10“ 4 

2.6 




600 

25 

0.21 

1.0 x 10 -4 

2.1 x 10 2 





100 

0.22 

1.1 x 10“ 4 

1.7 x 10 2 


...t 



100 

0.22 

1.1 x 10“ 4 

1.6 x 10 2 


...t 



100 

0.07 

2.8 x 10“ 5 

4.7 x 10 2 





200 

0.23 

1.1 x 10“ 4 

1.6 x 10 2 





400 

0.23 

1.1 x 10“ 4 

1.5 x 10 2 




1200 

100 

5.0 x 10“ 2 

1.8 x 10 -5 

8.3 x 10 3 


Core Galaxies, F = 

0.1 







10 6 


600 

25 

TI 

TI 

TI 




1200 

- 

1.5 x 10“ 2 

2.3 x 10“ 7 

1.7 x 10 s 


10 7 


600 

25 

0.19 

2.9 x 10 -5 

79 

0.6 




50 

TI 

TI 

TI 




1200 

25 

3.9 x 10“ 2 

1.4 x 10 -6 

2.9 x 10 4 





50 

4.1 x 10“ 2 

1.5 x 10 -6 

2.3 x 10 4 


10 8 


600 

25 

0.27 

5.7 x 10“ 5 

1.0 x 10 2 





50 

0.46 

1.5 x 10“ 4 

9.4 

0.2,0.6 



1200 

25 

8.5 x 10“ 2 

6.1 x 10 -6 

8.2 x 10 3 





50 

9.2 x 10“ 2 

7.0 x 10“ 6 

6.0 x 10 3 





100 

0.1 

8.8 x 10“ 6 

3.8 x 10 3 





200 

0.15 

1.7 x 10 -5 

8.0 x 10 2 



Break radius of stellar density profile. ^V; n f = 14as fixed in our numerical runs. T) Inflow rate in Eddington units, normalized 
to a stellar mass input parameter r/ = 0.2. { ' /!} Ratio of wind heating rate to radiative cooling rate at the stagnation radius. fCalculated 
with radiative cooling included, assuming mass loss parameter rj = 0.2. ^Calculated with radiative cooling and conductivity included, 
assuming mass loss parameter i] = 0.2 and conductivity saturation parameter <j> = 0.1. Solutions including radiative cooling were 
performed for cusp galaxies with v w = 300 km s _1 and core galaxies with v w = 600 km s —1 for tj = 0.02, 0.2, and 0.6. Values of ?/ 
resulting in thermally unstable solutions are marked in the final column. Solutions found to be thermally unstable for all r\ > 0.02 are 
denoted as TI. 


T(r), and radial velocity |v(r)|/c s , for the cusp (F = 0.8) 
solutions within our grid. As expected, the gas density in¬ 
creases towards the SMBH p oc r~ v with v ~ 1, i.e. shallower 
than the —3/2 power law for Bondi accretion. This power 
law behavior does not extend through all radii, however, 
as the gas density profile has a break coincident with the 
location of the break in the stellar light profile (rb = 100 
pc). The temperature profile is relatively flat at large radii, 
but increases as oc 1 /r k interior to the sphere of influence, 
where k < 1, somewhat shallower than expected for viri- 
alized gas within the black hole sphere of influence. The 
inwardly directed velocity increases towards the hole with a 


profile that is somewhat steeper than the local free-fall ve¬ 
locity v oc vs oc r -1,/2 . The flow near the stagnation radius 
is subsonic, but becomes supersonic at two critical points. 
The inner one at r ~ 0.1r s is artificially imposed for nu¬ 
merical stability, although we have verified that moving the 
inner boundary has a small effect on the solution properties 
near the stagnation radius. The outer one is located near 
the break radius, r ~ rb and is caused by the transition to 
a steeper stellar density profile exterior to the outer Nuker 
break radius. 

Figure 3 shows our calculation of the stagnation ra¬ 
dius r s /ri n f as a function of the wind heating parameter 
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Figure 2. Radial profiles of the CNM density (top), tempera¬ 
ture (middle), and velocity (bottom), calculated for a representa¬ 
tive sample of galaxies. Colors denote values of the effective wind 
heating rate, v w = 1200 km s -1 (blue), 600 km s —1 (orange), 
and 300 km s —1 (green). Line styles denote different black hole 
masses: M . = ! 0 6 Mg (dot-dashed), 10 7 Mq (solid), and 10 s Mg 
(dashed). Thin and thick lines denote cusp galaxies (T=0.8) and 
core galaxies (r=0.1), respectively. Squares mark the locations of 
the stagnation radius. 


£ = Q 1 + (v w /ao) 2 , with different colors showing different 
values of v w . Cusp and core galaxies are marked with square 
and triangles, respectively. Shown for comparison are our an¬ 
alytic results (eq. [A7]) with solid and dashed lines for cusp 
and core galaxies, respectively, calculated assuming v = 1 
and v = 0.6, respectively. 

Our analytic estimates accurately reproduce the numer¬ 
ical results in the high heating limit C, 2> 1 (v w 3> oo; 
r s < fmf). However, for low heating the stagnation radius 
diverges above the analytic estimate, approaching the stel¬ 
lar break radius r b ri n f. This divergence occurs approxi¬ 
mately when ( < Cc oc (rb/r S oi) 0 ' 5 ^ r-1) (eq. [16]). Physically 
this occurs because the heating rate is insufficient to unbind 
the gas from the stellar potential. Thus, for small values of 
the heating rate (small /) the location of the stagnation ra¬ 
dius will vary strongly with the break radius. This explains 
the behavior of the three vertically aligned green squares. 
These are three cusp (P = 0.8) galaxies with M . = 10 8 Mq 
and v w = 300 km s _1 but with different break radii (rb = 
25, 50, and 100 pc from top to bottom). This divergence of 
the stagnation radius to large radii occurs at a higher value 
of in core galaxies (T = 0.1), explaining the behavior of 
the two core galaxies shown in Fig. 3 as vertically aligned 
orange triangles. 

Figure 4 shows the inflow rate for a sample of our nu¬ 
merical solutions for different values of v w = 300, 600 and 
1200 km s" , and for both core (r=0.1) and cusp galaxies 
(r=0.8). Shown for comparison is our simple analytic esti¬ 
mate of M/A/edd from equation (18). For high wind veloci¬ 
ties (v w 3> cro) the stagnation radius lies well inside the black 
hole sphere of influence and our analytic estimate provides 
a good fit to the numerical results. However, for low wind 
velocities and/or high M . (large ao), the numerical accre¬ 
tion rate considerably exceeds the simple analytic estimate 
as the stagnation radius diverges to large radii (Fig. 3). 

Fig. 5 shows the ratio of wind heating to radiative 
cooling, gheat/|<jrad| (eq. [25] and surrounding discussion) 
as a function of radius for v w = 300,600, and 1200 km/s, 
M. = 10 7 Mq, and V = 0.8. Radiative cooling is calculated 
using the cooling function of Draine (2011) for solar metal- 
licity. For high heating rates of v w = 600 and 1200 km s _1 
cooling is unimportant across all radii, while for v w = 300 
km s _1 we see that qheat/|Qrad| and t coo i/ts can be less than 
unity, depending on the wind mass loss parameter r). Be¬ 
cause at the stagnation radius OWt/l^radl is within a factor 
of two of its minimum across the entire grid, the value of 
(f?heat /1 Qrad | )r s is a global diagnostic of thermal instability 
for cusp galaxies. For core galaxies, the minimum value for 
the ratio of wind heating to radiative cooling may be orders 
of magnitude less than the value at the stagnation radius. 
However, we find that the Q/C, c criterion (see eq. (16)) com¬ 
bined with the ratio of wind heating to radiative cooling at 
r s still gives a reasonable diagnostic of thermal stability for 
core galaxies. 

Also note that for v w = 600 and 1200 km s -1 we have 
?heat/|Qrad| ~ tcooi/iff near the stagnation radius, but these 
ratios diverge from each other at low heating (v w < 300 km 
s _1 ). This results because equations (19), (20) underesti¬ 
mate the true gas density in the case of subsonic flow (weak 
heating). 

Although most of our solutions neglect thermal conduc¬ 
tion and radiative cooling, these effects are explored explic- 
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Figure 3. Stagnation radius r s in units of the sphere of influ¬ 
ence radius ri n f (eq. [3]) for galaxies in our sample as a func¬ 
tion of the stellar wind heating parameter £ = y/l-\- ( v w /cro ) 2 . 
Green, orange, and blue symbols correspond to different values 
of v w = 300, 600, and 1200 km s —1 , respectively. Squares corre¬ 
spond to cusp galaxies (r = 0.8), while triangles correspond to 
cores (r = 0.1). Green circles correspond to cusp solutions which 
would be thermally unstable. The black curves correspond to the 
analytic prediction from equation (A7), with thick solid and dot- 
dashed curves calculated for parameters (r = 0.8, u ~ 1) and 
(r = 0.1, v ~ 0.6), respectively. The thin black solid line corre¬ 
sponds to the simplified analytic result for r s from equation (15) 
(recall that r[ n f ~ GM 0 

itly in a subset of simulations. Dagger symbols in Table 2 
correspond to solutions for which we turned on radiative 
cooling. For cases which are far from being thermal instabil¬ 
ity when cooling is neglected (e.g., M. = 10 8 M©, v w = 600 
km s _1 , r b = 100 pc, T] = 0.2, <j h eat/|<?rad| = 76), including 
radiative cooling has little effect on the key properties of the 
solution, such as the stagnation radius, mass inflow rate, and 
the cooling ratio gheat/|<?rad|- However, for solutions that are 
marginally thermally stable, including radiative losses acts 
to significantly decrease gheat/|<?rad|. We calculated solutions 
with radiative cooling for all cusp galaxies with v w = 300 
km s _1 and for all core galaxies with v w = 600 km s -1 , 
each for three different values of r] = 0.02,0.2,0.6, spanning 
the physical range of expected mass loss for continuous star 
formation (Fig. 7). Solutions found to be unstable for all 
values of r/ are marked in boldface as “TI” in Table 2, with 
the unstable values of rj provided in the final column. 

Including thermal conduction, by contrast, results in 
only order unity changes to our solutions for our fiducial 
value of the saturation parameter (f> = 0.1, consistent with 
our analytic expectations (eq. [23]). However, we found that 
for runs which included conductivity, the stagnation radius 
varied with the location of the inner grid boundary: as the 
boundary is moved inwards, more heat can be conducted 
outwards from deeper in the gravitational potential well 
of the black hole. We expect conductivity would begin to 
have a significant effect for an inner boundary < 0.01r s (the 
conductive run in Table 2 has an inner grid boundary at 



Figure 4. Inflow rate M/M e dd versus SMBH mass for galaxies in 
our sample, calculated for different values of the wind heating pa¬ 
rameter v w = 300 km s —1 (green), 600 km s -1 (orange), and 1200 
km s' 1 (blue). Squares correspond to cusp galaxies (T = 0.8), 
while triangles correspond to cores (r=0.1). The green circle cor¬ 
responds to a cusp solution which would be thermally unstable. 
Thin solid and dot-dashed curves correspond to our simple ana¬ 
lytic estimates of M/M e( jd (eq. [18]) for cusp and core galaxies, 
respectively. Thick curves correspond to the more accurate im¬ 
plicit analytic expression given by equation (A6). 



Figure 5. Ratio of the rates of heating to radiative cooling, 
9 heat /19rad I, as a function of radius (solid lines) for a M . = 
10 7 Mq cusp galaxy (T = 0.8). Dashed lines show the ratio of 
the cooling time-scale to the free-fall time-scale t coo i/tff. For high 
heating rates both ratios are approximately equal at the stagna¬ 
tion radius (squares). When 9hoat/9rad $ 10 (or, equivalently, 
icool/lff $ 10 near the stagnation radius), then the flow is sus¬ 
ceptible to thermal instabilities. 
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~ 0.03r s ). However, it is not clear physically if the magnetic 
field could remain coherent over so many decades in radius 
in the face of turbulence from stellar winds. 


5 HEATING SOURCES 


Key properties of the flow, such as the mass inflow rate 
and the likelihood of thermal instability, depend sensitively 
on the assumed heating rate oc qv^,. In this section we es¬ 
timate the heating rate, v w , taking into account contribu¬ 
tion from stellar winds, supernovae, millisecond pulsars and 
SMBH feedback. We make the simplifying assumption that 
all sources of energy injection are efficiently mixed into the 
bulk of CNM gas. In fact, slower stellar wind material may 
cool and form high density structures as in Cuadra et al. 
(2005). Additionally, some of the injected energy may leak 
away from the bulk of the CNM material through low den¬ 
sity holes, as described by Harper-Clark & Murray (2009) 
for stellar feedback and by Zubovas & Nayakshin (2014) for 
AGN feedback. 

We summarize the individual heating sources below, 
providing estimates of the value of the effective wind ve¬ 
locity 



for each, where e is volumetric heating rate. 


(32) 


5.1 Stellar winds 

Energy and mass input to the CNM by stellar winds is the 
sum of contributions from main sequence and post-main se¬ 
quence populations. At early times following star formation, 
energy input is dominated by the fast line-driven outflows 
from massive stars (e.g., Voss et al. 2009). At later times 
energy input is dominated by main sequence winds (e.g., 
Naiman et al. 2013). Mass input is also dominated by the 
massive stars for very young stellar populations, but for most 
stellar ages the slow AGB winds of evolved low-mass stars 
dominate the mass budget. 

In Appendix B we calculate the wind heating rate from 
stellar winds, v and the mass loss parameter, r/ (eq. [8]), 
as a function of age, t*, of a stellar population which is 
formed impulsively (Fig. Bl). At the earliest times (r+ < 10' 
yr), the wind heating rate exceeds 1000 km s -1 , while much 
later (r* ~ th) stellar wind heating is dominated by main 
sequence winds is much lower, v$, ~ 50 — 100 km s _1 . As 
will be shown in §5.5, for the case of quasi-continuous star 
formation representative of the average SFHs of low mass 
galaxies, the heating rate from stellar winds can also be 
significant, v „ ~ 1000 km s _1 . 

Stellar winds thus contribute a potentially important 
source of both energy and mass to the CNM. However, two 
additional uncertainties are (1) the efficiency with which 
massive stellar winds thermalize their energy and (2) heat¬ 
ing from core collapse supernovae, which is potentially com¬ 
parable to that provided by stellar winds. We neglect both 
effects, but expect they will act in different directions in 
changing the total heating rate. 
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5.2 Type la Supernovae 


Type la SNe represent a source of heating, which unlike core 
collapse SNe is present even in an evolved stellar population. 
If each SN la injects thermal energy Ei a into the interstel¬ 
lar medium, and the SN rate per stellar mass is given by 
.Ria, then the resulting volumetric heating rate of Eia Ria 
produces an effective heating parameter (eq. [32]) of 

la 

Vw =\j ---• ( 33 ) 

The thermal energy injected by each SN la, Ei a — ei a 10 51 
erg, depends on the efficiency ei a with which the initial blast 
wave energy is converted into bulk or turbulent motion in¬ 
stead of being lost to radiation. Thornton et al. (1998) es¬ 
timate a radiative efficiency ei a ~ 0.1, depending weakly on 
surrounding density, but Sharma et al. (2014) argues that 
eia can be considerably higher, ~ 0.4, if the SNe occur in a 
hot dilute medium, as may characterize the CNM. Hereafter 
we adopt ei a = 0.4 as fiducial. 

The SN la rate, Ri a , depends on the age of the stellar 
population, as it represent the convolution of the star forma¬ 
tion rate and the la delay time distribution (DTD) divided 
by the present stellar mass. In the limit of impulsive star 
formation, Ri a is the DTD evaluated at the time since the 
star formation episode. The observationally-inferred DTD 
(Fig. 1 of Maoz et al. 2012) has the approximate functional 
form 


Ria = 1.7 x 1(T 14 (r*/t h ) -112 Mg 1 yr- 1 (34) 


where r* is the time since star formation. 

From equations (33), (34) we thus estimate that 

Vw « 700(eia/0.4)°' 5 (r*/th)-°' 56 »?ao2 2 kms -1 

« 700(eia/0.4) o ' 5 (r*/th) 0 ' 09 kms -1 , (35) 

where the second line assumes r/ ~ 0.02(r*/th)- 1,3 for a 
single burst of star formation (e.g. , Ciotti et al. 1991). 

The high value of implies that la SNe represent an 
important source of CNM heating. However, SNe can only 
be approximated as supplying heating which is spatially and 
temporally homogeneous if the rate of SNe is rapid com¬ 
pared to the characteristic evolution time of the flow at the 
radius of interest (Shcherbakov et al. 2014). We define the 
“la radius” 

1/2 

~ 35M-g- 1 (r*/4)°' 56 pc (36) 



as the location exterior to which the time interval between 
subsequent supernovae ri a ~ (ARncRia) -1 ~ G/(ragRi a ) 
exceeds the local dynamical time-scale tg yn ~ r/ero, where 
we again adopt the la rate for an old stellar population and 
in the final equality we estimate ag ~ Q3a. using Al. — a 
(eq. [4]). 

Assuming that ao then by substituting 

(eq. 35) into equation (15) for the stagnation radius, we 
find that 


ria 

r s 


7%.j2 M .,8 1 ( £ ta/0-4)(r*/th) 056 core 

14 %02 M .“8' 1 ( e Ia/ 0 - 4 )( r *Ah) _0 ’ 56 CUSP 

7 M“g- 1 (ei a /0.4)(r*/4)°' 74 core 
14 A/-g' 1 (eia/0.4)(r*/th)°' 74 cusp 
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SN la can only be approximated as a steady heating source 
near the stagnation radius for extremely massive SMBHs 
with M, > 10 9 M 0 or for a very young stellar population 
with t+ < th- 

Even if SN la are rare near the stagnation radius, they 
may cap the inflow rate (and thus the SMBH accretion) rate 
by periodically blowing gas out of the nucleus of low mass 
galaxies. Between successive SNe, stars release a gaseous 
mass M g « yM^ri^/th interior to the la radius, which is 
locally gravitationally bound to the SMBH by an energy 
Ebind > Mgcrg. From the above definitions it follows that 


£la < (Vi *) 2 
Ebind 2<7 q 


(38) 


Hence, for low mass BHs with cto mJ/, SN la are 
capable of dynamically clearing out gas from radii ~ ri a > 
r a . Thus even when heating is sufficiently weak that the 
stagnation radius formally exceeds ri a , the SMBH accretion 
rate is still limited to a value 

Mia ^ yM. /n a \ 2 ~ r _ 

M e dd Meddth V r inf / 

J4.5 x 10 _ 4 M“g' 33 (r*/fh ) _0 ' 2 core 
( 2.2 x 10 _ 4 M“g' 84 (r*/th ) _0 ' 6 cusp , 1 ’ 

obtained substituting the la radius ri a (eq. [36]) for r a in the 
derivation leading to our estimate of M (eq. [18]). 

The deep gravitational potential wells of high mass 
galaxies prevent SN la from dynamically clearing out gas 
in these systems (<to 3> v^). Equation (39) nevertheless still 
represents a cap on the accretion rate in practice because the 
la heating rate (eq. [35]) is usually high enough to prevent 
the stagnation radius (calculated including the la heating) 
from substantially exceeding ri a . If the stagnation radius 
moves inwards from ri a , then decreased heating will force it 
outwards again. On the other hand, if the stagnation radius 
moves well outside of ri a , then the high level of la heating 
will force it inwards. 


5.3 Millisecond Pulsars 


Energy injection from the magnetic braking of millisecond 
pulsars (MSPs) is a potentially important heating source. If 
the number of MSPs per unit stellar mass is n map and each 
contributes on average a spin-down luminosity L a d, then the 
resulting heating per unit volume e « TsdU msp e map results 
in an effective heating rate (eq. [32]) of 


MSP 

Vw 


30 


/ Cmsp \l/2 ( L ad 
l 0.1 / \ 10 34 erg s _1 


1/2 

%)d 2 2 kms _1 , (40) 


where e map is the thermalization efficiency of the wind, nor¬ 
malized to a value < 0.1 based on that inferred by modeling 
the interstellar media of globular clusters (Naiman et al. 
2013). Our numerical estimate assumes a pulsar density 
Ums P ~3x 10 -4t> MSPs g _1 , calculated from the estimated 
~ 30, 000 MSPs in the Milky Way (Lorimer 2013) of stellar 
mass w 6 x 10 10 Mq. 

Based on the ATNF radio pulsar catalog (Manchester 
et al. 2005), we estimate the average spin-down luminosity 
of millisecond pulsars in the field to be L a< j ~ 10 34 erg s _ , 
resulting in u[)f SP < 30 km s _1 for y > 0.02. For higher 


spin-down luminosities, L a d — 10 35 ergs s -1 characteristic 
of some Fermi-detected pulsars, then the higher value of 
ttf SP < 300 km s " 1 makes MSP heating in principle impor¬ 
tant under the most optimistic assumptions e map = 1 and 
y = 0 . 02 . 

5.4 SMBH Feedback 

Feedback from accretion onto the SMBH represents an im¬ 
portant source of heating which, however, is also the most 
difficult to quantify (e.g., Brighenti & Mathews 2003, Di 
Matteo et al. 2005; Kurosawa & Proga 2009; Fabian 2012 
for a recent review). A key difference between AGN heating 
and the other sources discussed thus far is its dependence 
on the SMBH accretion rate M., which is itself a function 
of the heating rate (eq. [17]). 


5-4-1 Compton Heating 


There are two types of SMBH feedback: kinetic and radia¬ 
tive. Radiative feedback is potentially effective even in low 
luminosity AGN via Compton heating (e.g., Sazonov et al. 
2004, Ciotti et al. 2010), which provides a volumetric heat¬ 
ing rate (Gan et al. 2014) 

e = 4.1 x 10 - 35 n 2 £Tc ergcm -3 s _1 , (41) 

where £ = T/nr 2 is the ionization parameter and L is the 
SMBH luminosity with Compton temperature Tc ~ 10 9 K 
3 > T (e.g., Ho 1999, Eracleous et al. 2010). 

The importance of Compton heating can be estimated 
by assuming the SMBH radiates with a luminosity L = 
eMc 2 , where e is the radiative efficiency and where M is 
estimated from equation (17). Then using equations (15), 
(17), (20), (24) we calculate from equation (32) that the 
effective heating rate at the stagnation radius is given by 


, core 
,cusp, 


(42) 


where Tc ,9 = Tb/10 9 K and e _2 = e/0.01 ~ 1. We caveat 
that, unlike stellar wind heating, Compton heating depends 
on radius, scaling as vS(r) oc (n 2 £/p *) 1,/2 oc 
i.e. oc r -0 ' 6 and oc r -0 ’ 75 for core (r = 0 . 8 ; v ~ 1 ) and 
cusp (r = 0.1; v ~ 0.6) galaxies, respectively. Although our 
model’s assumption that the heating parameter be radially 
constant is not satisfied, this variation is sufficiently weak 
that it should not significantly alter our conclusions. 

If Compton heating acts alone, i.e. v w = v [/, then solv¬ 
ing equation (42) for 5 W yields 


, core 
,cusp. 


(43) 


Compton heating is thus significant in young stellar popula¬ 
tions with relatively high mass loss rates, e.g. vS > 300 km 
s _1 for y > 1 . 

The inflow rate corresponding to a state in which Comp¬ 
ton heating self-regulates the accretion flow is given by sub¬ 
stituting equation (43) into equation (18): 


M C _ f 2 x 10- 3 yS;Z 2 Tc° g 8 eI™ ,core 
M edd ~ \8 x 10 - 4 r 7 o.o 2 ° V°- 7 M.° 8 14 ,cusp. 


(44) 
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Sharma et al. (2007) estimate the value of the radiative 
efficiency of low luminosity AGN, e ra d, based on MHD shear¬ 
ing box simulations of collisionless plasmas. For Eddington 
ratios of relevance, their results (shown in their Fig. 6) are 
well approximated by 



opening angle 6§ =0.1 (characteristic of AGN jets) to prop¬ 
agate through a gaseous mass M g of radius r is estimated 
from Bromberg et al. (2011) to be 


fjet ~ 4000 yr 


h 


10 40 ergs" 1 


-1/3 


M s 


-T(- 

pc y V io 8 m 0 


1/3 


(49) 


Approximating M g ~ Mtg, the ratio of the jet escape time- 
scale to the dynamical time-scale tdyn ~ r/cr is given by 


where M, is the BH accretion rate. In general M, will be 
smaller than the inflow rate M calculated thus far by a factor 
a < 1 due to outflows from the accretion disc on small 
scales. Thus, the full efficiency relating the mass inflow rate 
to the radiative output is e = ae ra d, where we take a = 
0.1 following Li et al. (2013), who find that the fraction of 
the inflowing matter lost to outflows equals the Shakura- 
Sunyaev viscosity parameter of the disc. 

Substituting e into equation (44) and solving for the 
inflow rate results in 

M C f 4 x 10“ 2 r?o.o 2 T^° g 8 Mff 6 , core 

Medd ~ \9 x 10- 3 ^;g 2 Tc°' 7 M. 0 ; 8 14 ,cusp. 

Accretion may not be truly steady in cases where AGN 
heating dominates, due to the inherent delay between black 
hole feedback and the structure of the accretion flow on 
larger scales. Equation (46) may nevertheless represent a 
characteristic average value for the inflow rate if a quasi¬ 
equilibrium is achieved over many cycles. 

5-4-2 Kinetic Feedback 


£~7X10-X^)-'*, (50) 

independent of r. A jet with power sufficient to apprecia¬ 
bly heat the CNM on radial scales ~ r s (tj > 10~ 5 ) also 
necessarily has sufficient power to escape the nuclear region 
and propagate to much larger radii. Slower outflows from 
the accretion disc, instead of a collimated relativistic jet, 
provide a potentially more promising source of feedback in 
these systems. 

As in the case of Compton heating, kinetic heating 
could in principle ‘self-regulate’ the accretion flow insofar 
as a lower heating rate results in a higher accretion rate 
(eq. [48]), which in turn may create stronger kinetic feed¬ 
back. However, given the uncertainty in the efficiency of ki¬ 
netic heating, we hereafter neglect its effect and defer further 
discussion to §6.6. 

5.5 Combined Heating Rate 

The total external gas heating rate, 

V W = VM 2 + (<?SP)2 + („£)* + („• )2, (51) 


Kinetic feedback results from outflows of energy or momen¬ 
tum from close to the black hole in the form of a disc wind or 
jet, which deposits its energy as heat, e.g. via shocks or wave 
dissipation, over much larger radial scales (e.g. McNamara 
& Nulsen 2007; Novak et al. 2011; Gaspari et al. 2012a). 

Assume that the outflow power is proportional to the 
SMBH accretion rate, Tj = ejAf.c 2 = O.lejMc 2 , where 
ej < 1 is an outflow efficiency factor and we have again 
assumed a fraction a = 0.1 of the infall rate reaches the 
SMBH. Further assume that this energy is deposited as heat 
uniformly interior to a radius rheat and volume Vheat oc r 8 eat . 
The resulting volumetric heating rate e = O.lq Me 2 /Vheat 
near the stagnation radius results in a heating parameter 
given by 


v 


f 0.2t hej Me 2 \ 1/2 

\ Viieat 'HP* | r s ) 


(47) 


700 

V2MT 


kms 1 



where we have used the facts that M = r]M+\ re /th and 
p*(r a ) = M*(r s )(2 — F)/(47rr 8 ). If the bulk of the energy 
from kinetic feedback is released near the stagnation radius, 
then even a small heating efficiency q > 10 -4 is sufficient for 
u*, to exceed other sources of non-accretion powered heat¬ 
ing. However, if this energy is instead deposited over much 
larger physical scales comparable to the size of the galaxy, 
i.e. rheat > 10 kpc ~ 10 4 r s , then kinetic feedback is unim¬ 
portant, even for a powerful outflow with ej ~ 0.1. 

The time required for a jet of luminosity Tj and half 


includes contributions from stellar winds, supernovae, pul¬ 
sars, and radiative SMBH feedback. We implicitly assume 
the different sources of energy injection mix efficiently. In 
reality, this may not be the case as slower velocity sources 
cool and form high density structure due to high pressure 
from the environment (Cuadra et al. 2005). 

The strength of each heating source depends explicitly 
on the SMBH mass and the stellar population in the galactic 
nuclear region. The latter could best be described by a single 
starburst episode in the past, or by a more continuous SFH 
that itself varies systematically with the galaxy mass and 
hence M.. 


5.5.1 Single Starburst 

Figure 6 shows the contributions of heating sources as a 
function of time r* after a burst of star formation for black 
holes of mass M . = 10 6 M 0 (top) and M. = 10 8 M 0 (bot¬ 
tom). The heating and mass input parameters due to stellar 
winds, t)J,(r*) and ??(t*), are calculated as described in Ap¬ 
pendix B (Fig. Bl). The SN la and Compton heating rates, 
v [0 and v are calculated from equations (35) and (42), 
respectively, which also make use of r?(r*) as set by stellar 
winds. 6 SN la heating depends also on a convolution of the 
DTD (eq. [34]) and the SFH, which reduces to the DTD itself 

6 Both and ?;(/' (through r s /r i a ) depend on the total heating 
rate (eq. [51]), requiring us to simultaneously solve a series of 
implicit equations to determine each. 
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Figure 6. Sources contributing to the gas heating rate at time 
r* after a single burst of star formation. Top and bottom pan¬ 
els show black hole masses of M. = 10 6 Mq and M* = 10 8 Mq 
(both cusps with F = 0.8). Solid and dashed lines show the ranges 
of r* for which the accretion flow is thermally stable and unsta¬ 
ble, respectively, according to the ratio of <?heat /1 <?rad I near the 
stagnation radius (eq. [25]). Shown with horizontal gray lines are 
the stellar velocity dispersion for each SMBH mass, estimated as 
o-Q = a. (eq. [4]). 


for a single star burst. We account for the expected suppres¬ 
sion of SN la heating resulting from its non-steady nature 
by setting v ^ = 0 when the stagnation radius (calculated 
excluding la heating) exceeds the la radius ri a (eq. [36]). 

Figure 6 shows that stellar winds are the most impor¬ 
tant heating source at early times, r* < 10' years. Compton 
heating becomes more important at later times due to (1) 
the higher accretion rates that accompany the overall de¬ 
crease in all sources of heating, coupled with (2) the persis¬ 
tently high mass loss rates and gas densities associated with 
the still relatively young stellar population. Interestingly, for 
M. = 10 6 Mq, there are two consistent solutions for times 
10 8 yr < t* < 10 9 yr: one which has high la heating and 
low compton heating and one with low la heating and high 
compton heating. We choose the latter, but our model can¬ 
not truly distinguish between these scenarios. Also, neither 
scenario would produce a thermally stable flow. For high- 


Figure 7. Top Panel: Sources contributing to the total heat¬ 
ing rate of the CNM, v w (black): stellar wind heating (orange), 
la supernovae (green), millisecond pulsars (blue), and compton 
heating (pink). Each heating source varies with black hole mass 
M . calculated for average SFHs from Moster et al. (2013); see 
Appendix B for details. Bottom Panel: The ratio of the total 
heating rate (ijheat oc ) from the top panel to the radiative 
cooling rate (|<j ra d|) at the stagnation radius (red), parameter 
// characterizing stellar mass loss rate (eq. [8]) as a function of 
black hole mass M,, calculated for average star SFHs from Moster 
et al. 2013 (blue), the ratio of la radius to stagnation radius as 
a function of M, (green), and the ratio of C/Cc (purple), where 
£ c (eq. [16]) is the critical heating parameter ( = y/l + ( v w /ctq) 2 
below which outflows are impossible. All quantities are calculated 
for core (T = 0.1) galaxies. The results for cusp galaxies are qual¬ 
itatively similar. 

mass black holes (Fig. 6, bottom panel), SN la dominate at 
all times after 10 7 yr. Both Type la and Type II supernovae 
could be an important heating source for times r* < 40 Myr, 
albeit a non-steady one for M . < 10 8 Mg. For simplicity we 
neglect supernova heating for t* < 40, keeping in mind that 
it could raise the overall heating by a factor of a few. 

As the total heating rate declines with time, the flow 
inevitably becomes thermally unstable according to the cri¬ 
terion (f)he a t/[f/r a d|)r s 10 (eq. [27]), as shown by dashed 

lines in Fig. 6. Compton heating is neglected in calculat- 


MNRAS 000, 000-000 (0000) 

















Circumnuclear Media of Quiescent SMBHs 15 


ing thermal stability because—unlike local stellar feedback 
mechanisms—it is not clear that SMBH feedback is capable 
of stabilizing the flow given its inability to respond instan¬ 
taneously to local changes in gas properties. Thermally un¬ 
stable flow is present throughout the single starburst case 
except at very early times, r* < 10‘ yr and at late times 
r* > 10 9 yr in the high M. case. 

Finally, MSP heating is negligibly small for fiducial pa¬ 
rameters and hence is not shown, while kinetic feedback is 
neglected given its uncertain efficiency (§6.6). 


5.5.2 Continuous SFH 

A single burst of star formation does not describe the typ¬ 
ical star formation history of most galaxies. Qualitatively, 
smaller galaxies will on average experience more recent star 
formation, resulting in energetic young stellar winds and su¬ 
pernovae which dominate the gas heating budget. Massive 
galaxies, on the other hand, will on average possess older 
stellar populations, with their heating rates dominated by 
SN la (on large radial scales) and SMBH feedback. We esti¬ 
mate the average value of v w as a function of SMBH mass 
for each heating source by calculating its value using the av¬ 
erage cosmic star formation histories of Moster et al. (2013). 
Appendix B describes how the average SFH is used to deter¬ 
mine the stellar wind heating v„ and mass return parameter 
7 ] as a function of M. (Fig. 7; top panel). SFHs are also con¬ 
volved with the la DTD distribution to determine the SN 
la heating. 

Figure 7 shows v w (M.) from each heating source (stars, 
MSPs, SNe la, and black hole feedback) calculated for the 
average SFH of galaxies containing a given black hole mass. 
The younger stellar populations characterizing low mass 
galaxies with M. < 3 x 10 8 Mg are dominated by stellar 
winds, with v^ > 700 km s _1 . Only for M. > 3 x 10 8 Mg 
does the lack of young stellar populations significantly re¬ 
duce the role of stellar wind heating. For these massive 
galaxies, however, the la radius is sufficiently small that SN 
la contribute a comparable level of heating, v^ > 500 km 
s' 1 (SN la do not contribute to the heating in low mass 
galaxies because ri a > r s ; bottom panel). 

Strikingly, we find that galactic nuclei that experience 
the same average SFH as their host galaxies possess ther¬ 
mally stable flows across all black hole masses. An impor¬ 
tant caveat, however, is that v w is generally only a few 
times higher than the stellar velocity dispersion, i.e. C ~ Cc 
(eq. [16]). This implies that the true stagnation radius and 
the inflow rate could be larger than our analytic estimates, 
potentially resulting in thermal instability in a significant 
fraction of galaxies. The bottom panel of Figure 7 also shows 
C/Cc, where £ c is estimated using fits for rb, r;„f, and T de¬ 
rived from the Lauer et al. (2007) sample. In other words, the 
CNM of a significant fraction of massive ellipticals may be 
thermally unstable, not necessarily because they are receiv¬ 
ing lower stellar feedback than the average for their galaxy 
mass, but due to the high heating required for stability given 
the structure of their stellar potential. 

Realistic variations (“burstiness”) in the SFH can also 
produce lower stellar wind heating rates (Appendix B). Fig¬ 
ure B2 shows the heating rate as a function of average black 
hole mass for non-hducial cases in which the current (z = 0) 


star formation is suppressed by a factor i from its average 
z = 0 value, for a characteristic time-scale St *. For < 10' 
yr, we see that v^ is reduced by at most a factor of « 2 
from its average value, even for huge drops in the star for¬ 
mation rate (t ~ 10 -3 ). However, burstiness in the SFH 
over longer time-scales can suppress heating more signifi¬ 
cantly. When St* > 10 8 yr, v „ becomes very sensitive to u: 
at t = 0.1, v%, « 400 km s _1 in most galaxies, but for smaller 
l, « 200 km s -1 (an exception to this is in galactic nuclei 
with M. > 10 8 Mg, where heating is stabilized by SN la). 
Our use of average cosmic SFHs is appropriate provided lo¬ 
cal variations are either on short time-scales (<5f* < 10' yr), 
or limited in magnitude ((. > 0.1). If both of these conditions 
are violated, then the effective heating rates for all but the 
most massive galaxies (where la explosions dominate) will 
fall by a factor > 4 from our fiducial volumetric averages, 
calculated using SFHs from Moster et al. (2013). 

Another potential complication is the discreteness of 
local sources. We assume heating and mass injection are 
smooth, but energy and mass injection may be dominated 
by a handful of stars, particularly for small Al,. For example, 
for M. < 10 8 Mg, O stars dominate the energy injection in 
the average star formation histories of Moster et al. (2013), 
but the expected number of O stars inside the stagnation 
radius is 0.01 for Al. = 10 6 Mg and only exceeds 1 for M, > 
10' Mg, assuming circular stellar orbits 7 . 

6 IMPLICATIONS AND DISCUSSION 

6.1 Inflow and Black Hole Accretion Rates 

Figure 7 shows that the heating associated with the aver¬ 
age SFH of a galaxy is approximately constant with SMBH 
mass, except for the largest black holes with Al, > 5 x 
10 s Mg. For a fixed wind heating parameter, the Eddington 
ratio M/M e dd increases oc AfJ' 5 °' 8 ^ for core(cusp) galaxies, 
respectively (Fig. 4; eq. [18]). 

Figure 8 (top panel) shows with solid lines the inflow 
rate as a function of black hole mass for the average star 
formation heating, calculated from equation (18) using our 
results for the total wind heating (Fig. 4). This average in¬ 
flow rate increases from M ~ 10~ 6 — 10 _5 Mq yr -1 for low- 
mass black holes (M. ~ 10 6 Mg) to M ~ 10 _4 Mq yr -1 
for M. ~ 10 8 Mq. These fall below the maximum thermally 
stable inflow rate, Mti (yellow lines). 

For low mass black holes, blow out from SN la caps the 
SMBH accretion rate at a value Mi a ~ 10 -6 — 10~ 5 Mg yr -1 
(green lines; eq. [39]), which is however not low enough to 
prevent otherwise thermally-unstable flow at smaller radii. 
For higher M. the stagnation radius resides close to the la 
radius so la heating contributes to the steady gas heating 
rate, even if the energy released by a single la is insufficient 
to unbind gas from the stellar bulge (eq. [38]). This explains 
why the black and teal lines meet at high Al,. 

Shown also for comparison in Figure 8 (top panel) is the 
inflow rate, M, onto SgrA* calculated by Quataert (2004) 
and by Cuadra et al. (2008). Also shown is the range of M 

7 If the stars providing the heating reside on elliptical orbits, a 
greater number will contribute at any time due to the portion of 
their orbital phase spent at small radii. 


MNRAS 000, 000-000 (0000) 


16 Generozov, Stone, & Metzger 




Figure 8 . Top panel: Gas inflow rate A/as a function 
of black hole mass M, , shown for core (T = 0.1, dot-dashed) 
and cusp galaxies (T = 0.8, solid). Black lines show the inflow 
rate calculated from equation (18) using the heating rate pro¬ 
vided at z = 0 by the average star formation histories for each 
galaxy mass (Fig. 7). Teal lines show the maximum accretion 
set by la supernovae blow-out or heating, Mi a (39). Red lines 
show the inflow rate obtained if Compton heating acts alone, 
Me (eq. [46]). Yellow lines show the maximum inflow rate for a 
thermally stable flow near the stagnation radius, Mti (eq. [29]). 
Also shown are the inferred inflow rates of SgrA* from Quataert 
(2004) and Cuadra et al. (2008) (gray circles) and for NGC3115 
from Shcherbakov et al. (2014) (gray diamonds). Bottom panel: 
Growth times t grow = M,/M, in units of the Hubble time / ^ 
for each of the accretion rates shown in the top panel, where 
M, = aM, where a = 0.1, to account for the fraction of inflow¬ 
ing mass lost to disc outflows on small scales. 

for the low-luminosity AGN NGC3115 derived through de¬ 
tailed modeling by Shcherbakov et al. (2014), who find a 
range of inflow rates depending on the assumed model for 
thermal conductivity. The inflow rate SgrA* estimated from 
Quataert (2004) 8 exceeds our estimate of M due to the aver¬ 
age SFH for the same M. by at least an order of magnitude. 
This is because the mass source term in Quataert (2004) ex- 

8 There is a typo in Quataert (2004): the inflow rate corresponds 
total stellar mass loss rate of 5 X 10 _4 Mq yr —1 not 10" 3 Mq 
yr — , as pointed out by Cuadra et al. (2006). 


ceeds ours by a factor of ~ 50 (for cusp galaxies). The star 
formation in the central parsec of the Milky Way cannot be 
described as steady-state: feedback is instead dominated by 
the stellar winds from the ring of young massive stars of 
mass ~ 10 4 Mq and estimated age ~ 10 Myr (e.g., Schodel 
et al. 2007). Such a star formation history is better described 
by our impulsive (starburst) scenario, for which the value of 
r/ ~ 30 at t* = 10 7 yr (Fig. Bl) is a factor of ~ 100 times 
higher than the value of r/ ~ 0.4 predicted in the average 
SFH case (Fig. 7, bottom panel). Note that the maximum 
thermally stable accretion rate is dependent on the r/ param¬ 
eter and would be higher in the case of an impulsive burst 
of star formation. 

However, simulations by Cuadra et al. (2006) find that 
the inflow rate for Sgr A* is reduced by up to an order 
of magnitude compared to the result in Quataert (2004) 
when the discreteness and disc geometry of stellar sources 
is taken into account, and recent results from Yusef-Zadeh 
et al. (2015) suggest that the inflow rate could be further 
reduced to M / Mem ~ 3 x 10 -6 once clumping of stellar 
winds is taken into account. 


6.2 Nuclear X-ray Luminosities 

The unresolved core X-ray luminosities of nearby galactic 
nuclei provide a powerful diagnostic of the SMBH accretion 
rates and how they vary with SMBH mass and other galaxy 
properties (e.g., Ho 2008 and references therein). Figure 9 
shows our predictions for {Lx) = ex Me 2 as a function of 
black hole mass, where ex = cue ra d£boi, cx = 0.1 accounts for 
the fraction of the inflowing matter loss to outflows from the 
accretion disc, e ra d is the radiative efficiency of the accreted 
matter, and tboi = 0.1 is the assumed bolometric correction 
into the measured X-ray band for low luminosity AGN (Ho 
2008). In all cases the value of {Lx) is calculated using the 
accretion rate from the top line of equation (17) with F = 0.7 
for M. < 4 x 1O 7 M 0 , and T = -O.31og 10 ( 4x ™’ Mr . ) +°- 7 

for M m > 4 x 10' M 0 . This functional form is designed to 
approximately reproduce the behavior of F(M.) in the Lauer 
et al. (2007) sample. 

The left panel of Figure 9 is calculated assuming a con¬ 
stant low value for the radiative efficiency of ex = ICG 4 , 
typical of those estimated for low luminosity AGN (e.g., 
Ho 2009). The right panel luminosities are calculated in¬ 
stead assuming an efficiency of ex = eboiGad that depends 
on the Eddington ratio as predicted by MHD shearing box 
simulations by Sharma et al. (2007, see eq. [45] and sur¬ 
rounding discussion). Shown for comparison are the X-ray 
measurements (black stars) and upper limits (gray trian¬ 
gles) from the sample of early-type galaxies compiled by 
Miller et al. (2015, cf. Gallo et al. 2010). A black line shows 
the best power-law fit to the X-ray luminosity from Miller 
et al. (2015), given by {Lx/L b m) ~ M , with a = —0.2 (see 
also Zhang et al. 2009; Pellegrini 2010; Gallo et al. 2010). 
Also shown are the maximum accretion rates, respectively, 
for thermally stable accretion (eq. [29]; yellow), as set by 
SN la blow-out/heating (eq. [39]; teal), and as allowed by 
Compton heating feedback (eq. [46]; pink). 

If the nuclei of elliptical galaxies are heated as expected 
for the average SFH of galaxies with similar mass, then to 
first order we predict that {Lx/L e dd) should be an increas- 
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Figure 9. Average nuclear X-ray luminosity, (Lx) = ex Me 2 , as a function of SMBH mass. Green lines show our prediction in the 
case of heating due to the average SFH for galaxies corresponding to each black hole mass (Fig. 8), calculated using T = 0.7 for 
M,<4x 10 7 M o , and T = -O.31og 10 (M./4 x 10 7 M©)+0.7 for M, >4x 10 7 M Q , as derived from the Lauer et al. (2007) sample. 
Shown for comparison are the measurements ( black stars ) and upper limits ( gray triangles) for Lx/Ledd values, from the Miller et al. 
(2015) sample of early-type galaxies (a black line shows the best power-law fit (Lx/L e dd) oc M“, with a = —0.2). The left panel is 
calculated assuming a constant radiative efficiency ex = 10 — 4 (yi K) i/(). 1), while the right hand panel assumes ex = »CbolCrad> where e ra d 
is the radiative efficiency of low luminosity accretion discs calculated by Sharma et al. (2007) (eq. [45]) and a = 0.1 is the fraction of the 
mass inflowing on large scales that reaches the SMBH. Shown for comparison are the X-ray luminosities calculated for the maximum 
thermally-stable accretion rate ( dashed orange line ; eq. [29]); the SN la-regulated accretion rate ( dashed teal line-, eq. [39]), and the 
Compton heating-regulated accretion rate ( dashed pink line-, eq. [44] ) 10 , again all calculated for the average SFH corresponding to each 
SMBH mass. 


ing function of BH mass. This result is in tension with the 
observed down-sizing trend: the average Eddington ratio is 
seen to decrease (albeit weakly) with A1., especially in the 
model where ex increases with the Eddington ratio. How¬ 
ever, one must keep in mind the enormous uncertainty in 
calculating the luminosity of the accretion flow close to the 
black hole in a single waveband based on the feeding rate 
on larger scales. Also, if disc outflows indeed carry away 
most of the infalling mass before it reaches X-ray producing 
radii (e.g. Blandford & Begelman 1999; Li et al. 2013), then 
the angular momentum of the infalling gas might also influ¬ 
ence the X-ray luminosity indirectly through the inflow effi¬ 
ciency a, in a way that could depend systematically on the 
stellar population and hence M,. However, in general the 
efficiency with which inflowing matter reaches the SMBH 
(and hence X-ray luminosity) would naively be expected to 
decrease with decreasing At, due to the increasing angular 
momentum of low mass galaxies, exacerbating the tension 
between our findings and the observed downsizing trend. 

Perhaps a more readily addressable test is whether 
a steady-state accretion picture developed in this work 
can account for the the large scatter, typically of 2 — 3 
orders of magnitude, in Lx/L e dd at fixed Al.. Scatter 
could result from the strong sensitivity of the accretion 
rate to the stellar wind velocity and mass loss parameter, 
M/M e dd oc u^ 3 ' 8(_2,4) for core(cusp) galaxies, respectively 
(Fig. 4; eq. [18]). When combined with the significant depen¬ 
dence of v w on stochastic intermittency in the star formation 
history (Appendix B2), this can lead to order of magnitude 
differences in At. However, we note that v w is not expected 


to vary by more than a factor of a few for a thermally-stable 
flow, limiting the allowed variation of M/M e dd- The discrete 
nature of stellar wind sources and their motions could result 
in an additional order of magnitude variation At (Cuadra 
et al. 2008). 

Variations in Lx / L e dd could also result from differences 
in angular momentum of the infalling gas from galaxy to 
galaxy, resulting in differences in the fraction of the gas 
lost to outflows. Also potentially contributing is the order 
of magnitude difference, at fixed v w and At., between M for 
core and cusp galaxies (Fig. 8). Differences in M are aug¬ 
mented by the theoretical expectation that Lx oc At 2 for 
radiatively inefficient flows. 


6.3 Steady Accretion versus outbursts 

It is also possible, and indeed likely, that many of the Miller 
et al. (2015) sample of galactic nuclei are not accreting in 
steady state. This is supported by the fact that many of the 
X-ray luminosities in Figure 9 lie above the predictions for 
stable accretion, i.e., Lx > exALric 2 (yellow line), at least 
for our choice of radiative efficiency. 

Within our steady state solutions, gas outside the stag¬ 
nation radius r s is blown out in wind. However, for very 
massive galaxies (At. > 10 8 Mg), the energy injected into 
the gas is not enough to truly eject it from the potential 
well of the galaxy, even if the gas is blown out of the stellar 
bulge (e.g. even if the C > Cc criterion is satisfied-see equa¬ 
tion (16) and surrounding discussion). This gas is likely to 
build up on the outskirts of the galaxy until a thermal in- 
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stability develops, causing the CNM to undergo cyclic oscil¬ 
lations of rapid cooling and high accretion rates, followed by 
quiescent periods once gas has been consumed and/or AGN 
feedback becomes effective. Such cycles are seen in numeri¬ 
cal simulations of elliptical galaxies on large scales and long 
time-scales (Ciotti et al. 2010). Evidence for such periodic 
outbursts includes observations showing that a fraction of 
early-type galaxies in the local Universe have undergone re¬ 
cent (< 1 — 2 Gyr old) star formation episodes (Donas et al. 
2007). 

If the actual radiative efficiency is lower than our fidu¬ 
cial assumption, then an even larger fraction of the X-ray 
detections would lie above the predictions for stable accre¬ 
tion. For example, we assume that a fraction a = 0.1 of the 
mass inflow rate reaches the SMBH, but this fraction could 
be considerably smaller. In such a case the X-ray detected 
galaxies will experience less heating than we calculate for the 
average SFH, and may not be stably accreting at all. This is 
plausible given that the Miller et al. (2015) sample includes 
only early-type galaxies with older stellar populations than 
the average for their stellar mass. 

Figure 8 (top panel) also shows that Me > Mri across 
all black hole masses. This indicates that Compton heating 
is not sufficient to produce a thermally stable accretion rate, 
at least for conditions corresponding to the average SFH. 

Similar cyclic AGN activity could occur on the smaller 
radial scale of the sphere of influence (e.g. Yuan & Li 2011, 
Cuadra et al. 2015). In this case the large scatter in the X- 
ray luminosities shown in Fig 9 could result from the wide 
range of inflow rates experienced over the course of a cyclic 
episode between instabilities and periodic nuclear outbursts. 

If the galaxies with measured Lx at low M, in Fig¬ 
ure 9 are indeed in a state of thermally-unstable outburst, 
then some of the galaxies with X-ray upper limits could be 
contributing to a separate population of thermally-stable, 
steadily accreting nuclei. Such a population could include 
SgrA*, which has a much lower X-ray luminosity for its 
SMBH mass than predicted by the Miller et al. (2015) trend. 
A potentially bimodal population of steady and outbursting 
galactic nuclei calls into question the practice of using simple 
extrapolations of power-law fits to the Lx(M.) relationship 
to low M . to constrain the occupation fraction of SMBHs 
in the nuclei of low mass galaxies. 

6.4 SMBH Growth Times 

Low-mass AGN in the local Universe with M. < 3 x 10 7 Mq 
are observed to be growing on a time-scale comparable to 
the age of the Universe, while the most massive SMBHs with 
M. > 10 9 Mq possess local growth times which are more 
than 2 orders of magnitude longer (Heckman et al. 2004; 
Kauffmann & Heckman 2009). 

The bottom panel of Figure 8 shows the SMBH growth 
time, tgrow = M./M for each accretion rate shown in the 
top panel, recalling that M. = aM , where a = 0.1. For 
SMBHs which accrete steadily at the rate set by stellar wind 
heating due to the average star formation history of their 
host galaxies, we see that f grow exceeds the Hubble time 
by 2—4 orders of magnitude across all M,. Steady accretion 
therefore cannot explain the growth of low mass black holes, 
a fact which is not surprising given that approximately half 
of this growth occurs in AGN radiating within 10 per cent 


of their Eddington rate (Heckman et al. 2004). Such high 
accretion rates likely instead require a source of gas exter¬ 
nal to the nuclear region, triggered either by galaxy mergers 
associated with the hierarchical growth of structure or ther¬ 
mal instabilities on larger, galactic scales (e.g. Ciotti et al. 
2010, Voit et al. 2015). 

That the growth time associated with thermally- 
unstable accretion (yellow line in Fig. 8 ) exceeds the Hubble 
time across all SMBH masses highlights the fact that signif¬ 
icant black hole growth in the local Universe cannot result 
from thermally-stable steady accretion of gas lost from the 
surrounding stellar population studied in this paper. Gas 
blow-out by SN la cannot alone prevent the growth of low- 
mass black holes, as indicated by the low growth times <C th 
allowed by la heating (green lines), although la heating 
could play in principle a role in capping the growth rate 
of > 1O 8 M 0 black holes, again depending on the efficiency 
of la heating. 


6.5 TDE Jets 

Our results also have implications for the environments 
encountered by relativistic jets from TDEs. For low-mass 
SMBHs with M. < 10 'Mq, such as that responsible for 
powering the transient Swift J1644+57 (Bloom et al. 2011), 
we predict gas densities of a n ~few cm -8 on radial scales of 
0.3 pc for wind heating rates within the physically expected 
parameters u w ~ 500—1000 km s -1 and 77 = 0.4 (see Fig. 2). 
This is comparable to the density ~ 0.3-10 cm -3 obtained 
by modeling the radio afterglow of Swift J1644+57(Berger 
et al. 2012, Metzger et al. 2012). However, there is consid¬ 
erable uncertainty in the afterglow modeling. For example, 
Mimica et al. (2015) find a much higher density of 60 cm ~ 3 
at 0.3 pc. 

Berger et al. (2012) infer a flattening of the gas density 
profile of the host of J1644+57 on radial scales r > 0.4 pc 
(their Fig. 6 ) which looks qualitatively similar to the shape 
of the density profile we predict for core galaxies (Fig. 2, 
solid line). To obtain such a flattening on the inferred radial 
scales would require a black hole of M, ~ 10' Mq. However, 
this is inconsistent with the break in the gas density seen at 
~ 0.6 pc by Berger et al. (2012), which should correspond to 
the break radius of the stellar density profile, thus implying 
an unrealistically small SMBH mass of M. ~ 100 Mq ac¬ 
cording to the mean r^ — M, relationship (Lauer et al. 2007). 
Alternatively, the observed density break in J1644+57 could 
result from the outer edge of a sub-parsec nuclear star clus¬ 
ter (e.g., Carson et al. 2015). The high stellar densities of 
such a compact star cluster could greatly enhance the TDE 
rate (e.g. Stone & Metzger 2014), possibly making this as¬ 
sociation uncoincidental. On the other hand, we note that 
some models for J1644+57 (e.g., Tchekhovskoy et al. 2014) 
favor a much smaller black hole (10 5 — 10 6 Mq) than we 
estimate would be required to produce the observed inner 
break. 
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6.6 Regulation by SMBH Kinetic Feedback 

Our analysis of SMBH feedback has focused on Compton 
heating instead of kinetic feedback, since the effects of ra¬ 
diative feedback are relatively straightforward to calculate 
from the properties of the accretion flow. However, it is pos¬ 
sible kinetic feedback from SMBH outflows could play an 
equal or greater role in regulating accretion, even in low 
luminosity AGN. 

Our simple parameterization of kinetic feedback 
(eq. [48]) assumes a uniform volumetric heating and pre¬ 
dicts an effective wind velocity which increases with SMBH 
mass as vj, oc Mf 88 for core and cusp galaxies, respectively. 
Coupled with the dependence of the SMBH accretion rate on 
the wind heating, M/M e dd oc Mf 7bl ' 0A8 ' > v^, 3 ' 81, 2,4) (Fig. 4; 
eq. [18]), a dominant source of kinetic feedback of the form 
we have adopted leads to an Eddington ratio dependence on 
SMBH mass of M/M edd oc M~°- 7{ ~ 0A) . 

This prediction is more in line with the observed de¬ 
pendence of (Lx/L e dd) in elliptical galaxies (Fig. 9, green 
line), suggesting that accretion regulation by kinetic feed¬ 
back could play a role in determining the X-ray luminosities 
of elliptical galaxies. However, our assumption that kinetic 
outflows heat the gas uniformly in volume is rather arbitrary 
and would need to be more rigorously justified by numerical 
studies of how jets or disc outflows couple energy to their 
gaseous environment. As already discussed, it is furthermore 
unclear whether a steady-state accretion model is at all rel¬ 
evant in the case when SMBH feedback dominates due to 
the time delay between the small scale accretion flow and 
the thermodynamic response of the CNM on larger scales. 


7 CONCLUSIONS 

We have calculated steady-state models for the hot gaseous 
circumnuclear media of quiescent galaxies, under the as¬ 
sumption that gas is supplied exclusively by stellar wind 
mass loss and heated by shocked stellar winds, supernovae 
and black hole feedback. We numerically compute solutions 
for a range of different black hole masses, heating rates, and 
observationally-motivated stellar density profiles. Then we 
use our numerical results (Table 2) to verify and calibrate 
analytic relationships (Appendix A). We use the latter to 
explore systematically how the SMBH accretion rate varies 
with black hole mass and the galaxy’s SFH. Our results for 
M(M.) are compared with observed trends of the nuclear 
X-ray luminosities of quiescent SMBHs and low luminosity 
AGN. Our conclusions are summarized as follows. 

(i) A stagnation radius, r a , divides the nuclear gas be¬ 
tween an accretion flow and an outgoing wind. In steady- 
state the gas inflow rate towards the black hole is propor¬ 
tional to the stellar mass enclosed inside of r s . In the limit of 
strong heating, the stagnation radius resides interior to the 
SMBH influence radius and coincides with the Bondi radius 
(eq. [22]). In the limit of weak heating (£ < eq. [16]), the 
stagnation radius moves to large radii, near or exceeding the 

10 These equation correspond to the specific cases of T = 0.1 and 
T = 0.8, but it is straightforward to generalize these for arbitrary 
T as we do for this figure. 


stellar break radius r d ri n f, greatly increasing the density 
of gas on smaller radial scales. 

(ii) In the vicinity of the stagnation radius, including the 
effects of heat transport by electron conduction results in 
at most order unity changes to the key properties of the 
flow (e.g. stagnation radius, mass accretion rate, and ther¬ 
mal stability) for causal values of the conduction saturation 
parameter (j> < 0.1 (eq. [23]). However, in principle heat 
conduction of heat from the inner accretion flow can affect 
the solution properties on much larger scales (Johnson & 
Quataert 2007). For example, Tanaka & Menou (2006) find 
that conductivity can contribute to driving bipolar outflows. 

Angular momentum of the gas will become important on 
small scales, where gas stopped by a centrifugal barrier and 
unable to cool could drive outflows near the equatorial and 
polar regions of the flow. 

(iii) Radiative cooling has a more pronounced influence 
on the flow structure when the radiative cooling rate ex¬ 
ceeds the gas heating rate near the stagnation radius, where 
the gas is in nearly hydrostatic equilibrium (Fig. 1). This 
condition is approximately equivalent to the i coo i < lOtfi cri¬ 
terion for thermal instability advocated by e.g. Sharma et al. 
(2012), Gaspari et al. (2012b) and Li & Bryan (2014a). The 
transition in the flow properties that occurs as heating is re¬ 
duced may represent a true “thermal instability” (hot ISM 
condensing into a cooler clouds), or it may simply represent 
an abrupt transition from a steady inflow-outflow solution 
to a global cooling flow. We leave distinguishing between 
these possibilities to future work. 

(iv) The location of the stagnation radius, and hence 
the inflow rate, is a sensitive function of the gas heat¬ 
ing rate oc v%,. Quantitatively, r s oc v^ 2 , implying that 
M oc Vw 3 ' 8i '~~ A ' > for fiducial core (cusp) galaxies. However, 
M can increase even more rapidly with decreasing v w for 
two reasons: (1) when v w becomes comparable to the stellar 
velocity dispersion (£ < £ c ; eq. [16]), then gas remains bound 
to the stellar bulge and cannot produce an outflow interior 
to the stellar break radius (2) For low v w , gas near the stag¬ 
nation radius becomes thermally unstable, which based on 
the results of previous numerical studies (e.g.,Ciotti et al. 
2010 ) is instead likely to result in a burst of high accretion 
(§6.3). 

(v) Stellar wind heating, supernovae, and AGN feedback 
depend explicitly on the SMBH mass as well as the SFH 
in the galactic nucleus. A young starburst of age < 10 Myr 
produces mass and energy injection dominated by young 
stellar winds (v w ~ 1000 km s“ ) and supernovae, while 
for exclusively old stellar populations mass input is dom¬ 
inated by AGB winds and heating is dominated by AGN 
feedback and supernovae. Ongoing, continuous star forma¬ 
tion presents a hybrid situation, with energy input usually 
dominated by young stars and mass input dominated by old 
stars (Appendix B). Galactic nuclei with M, < 10 8 Mq that 
are heated according to their average SFHs (as derived from 
cosmological simulations) receive a stellar wind heating of 
v w ~ 700 km s" . This heating can be suppressed if the star 
formation is sufficiently intermittent (Fig. B2). 

(vi) Type la SNe only provide a continuous heating source 
exterior to the la radius ri a (eq. [36]) where the time between 
subsequent SNe exceeds the dynamical time-scale. Non-la 
heating due to the average SFH results in r B < ri a , except 
for the most massive galaxies (Fig. 7). However, if r s does 
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approach ri a , then SN la heating is usually large enough to 
prevent r a from greatly exceeding ri a . SNe la thus regulate 
the SMBH accretion rate below the value A/i a (eq. [39]), 
except possibly in high mass elliptical galaxies with large 
break radii and C( v i?) > Cc- 

(vii) Unlike stellar heating, heating from SMBH feedback 
depends on the accretion rate, which itself depends on the 
heating. Compton heating is generally unimportant com¬ 
pared to stellar wind heating for the average SFH, but it can 
be significant in the case of an impulsive starburst(Fig. 6). 
However, SMBH feedback may not be capable of truly stabi¬ 
lizing the flow given its inability to respond instantaneously 
to local changes in the properties of the gas. Kinetic feedback 
could also be important in determining M(M ,) when stellar 
heating is weak (§6.6) but is more challenging to quantify. 

(viii) Stellar wind heating from the average SFH may be 
sufficient to permit thermally-stable, steady-state accretion, 
depending on the accretion efficiency of the SMBH. How¬ 
ever, the resulting average Eddington-normalized accretion 
rate is predicted to increase with M. (Fig. 9), in tension 
with the (also weak) downsizing trend of measured X-ray 
luminosities in early-type galactic nuclei (e.g., Miller et al. 
2015). Thermally stable accretion models can reproduce the 
observed scatter in nuclear X-ray luminosities at fixed M. 
(two to four orders of magnitude) due to the combination 
of (1) differences in the stellar wind heating rate due to 
stochastic variation in SFH between galaxies, or within one 
galaxy due to burstiness in the SFH, as in Fig. B2; (2) 
variations in the amount of inflowing mass which is lost to 
small scale outflows from the SMBH accretion disc, likely 
due to variations in the angular momentum of the accreting 
gas (cf. Pellegrini 2010). 

(ix) However, for lower X-ray radiative efficiencies, the ac¬ 
cretion rates of early-type galaxies are above the maximum 
value for thermally stable flow, Mti (Fig. 9). This implies 
that current X-ray detections could instead be comprised 
mostly of nuclei undergoing outburst due to thermal insta¬ 
bilities. In such a case, there may exist a separate (usually 
undetected) branch of low-Lx nuclei accreting stably. The 
existence of such a bimodal population would call into ques¬ 
tion constraints on the SMBH occupation fraction in low 
mass galaxies (Miller et al. 2015) derived by extrapolating 
the Lx(M .) relationship to small M.. 

(x) Low mass black holes grow in the low redshift Uni¬ 
verse over time-scales comparable to the Hubble time (Heck¬ 
man et al. 2004). The accretion rates so required are too high 
to be consistent with either the value predicted for the aver¬ 
age SFH history or the maximum allowed for steady-state, 
thermally stable accretion (Fig. 8). Perhaps unsurprisingly, 
low -M, growth and AGN activity must instead be driven by 
a supply of gas external to the nucleus, such as galaxy merg¬ 
ers or thermal instabilities on larger, galactic scales (e.g. Voit 
et al. 2015). 

We conclude by drawing attention to a few limitations 
of our model. First, we have neglected any large scale inflow 
of gas to the nucleus by assuming the only source of gas 
feeding the black hole is stellar wind mass loss from the local 
stellar population. Our estimates represent a lower bound 
on the time averaged gas density of the CNM. Although our 
model can, under certain assumptions, describe quiescent 
galactic nuclei, it cannot account for AGN. 


Predictions for the SMBH accretion rate are strongly 
tied to the age and radial distribution of the stellar popu¬ 
lation within the sphere of influence. Our model therefore 
cannot make definitive predictions for the accretion mode of 
individual galaxies without knowledge of their specific SFH. 
We adopt halo-averaged star formation rates measured using 
multi-epoch abundance matching models by Moster et al. 
(2013), and neglect any spatial variation in the star for¬ 
mation rate for any given galaxy. However, galaxies of mass 
> 7x 10 10 Mq assemble their stars from the inside out (Perez 
et al. 2013), resulting in galactic nuclei which may be sys¬ 
tematically older than assumed in our model and thus more 
prone to thermal instability. 

Finally, we do not account for the effects of non- 
spherical geometry or discreteness of stellar wind sources 
(Cuadra et al. 2006, 2008). We assume all of the mass and 
energy sources are efficiently mixed, and do not take into 
account the possibility of slower wind material condensing 
into high density structures (as described in Cuadra et al. 
2005) or some of the energy injection leaking away in 3-D 
(Harper-Clark & Murray 2009; Zubovas & Nayakshin 2014). 
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APPENDIX A: DERIVATION OF STAGNATION 
RADIUS 


Here we derive an analytic expression for the time- 
independent stagnation radius, r a . In steady-state, the re¬ 
quirement that the net heating rate on the right hand side 
of the entropy equation (eq. [7]) must equal zero at the stag¬ 
nation radius (where v = 0 ) fixes the temperature at r 3 : 


7 P _ fv|r a khT\r s _ 7 1 f+|r s 

7-1 P r B 2 frnip 7 2 


(Al) 


Combining (eq. [7]) with the first law of thermodynamics, 


T ds 
dr 


1 

d(p/p) 

V p 

H-— 

r s P 

r s b r 

Q (V_ _L ^- ^P) 

2 T 2 7-I p/ 

r s 7-1 

dr 

ra pv 


(A2) 


where v = — d Inp/dlnrj . The right side can be evaluated 
using L’Hopital’s rule, yielding 
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where we have used the definition v%, = v^, + GA f' + 
ctq (eq. [9]). Substituting this expression back into equation 
(A2) and using equation (Al) gives 
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Combining this with condition of hydrostatic equilibrium 
(eq. [ 6 ]) at the stagnation point, 
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results in the following expression for the stagnation radius: 
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where we have assumed 7 = 5/3. This relationship can be 
reparameterized as 
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where ri n f = 3GM,/oo and we have defined £ = 

\/T+ {v w /a 0 ) 2 . Because M* | ra = M.(r s /n n i) 2 ~ r , in general 
equation (A7) must be solved implicitly for r s /r m {. However, 
when M*| ra << M. or equivalently r s 4C ri n f, equation (A7) 
simplifies considerably: 
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APPENDIX B: ANALYTIC MODEL FOR 
DEPENDENCE OF WIND HEATING v* ON 
STELLAR AGE 

B1 Single Burst 


In the case of a single impulsive burst of star formation, the 
mass and energy injection rate per unit stellar mass at time 
t after the burst are given, respectively, by 
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where the f’s represent the efficiency with which each 
source of mass/energy injection is thermalized (we take all 
of the f’s to be 1). In the top line, the first term corresponds 
to mass injection due winds from post-main sequence stars. 
The second term corresponds to mass injection due to stellar 
winds from massive stars. For the latter we use population 
synthesis calculations by Voss et al. (2009). From the bottom 
panel of their Figure 7, the stellar wind mass loss rate per 
massive star can be approximated to within a factor of a few 
by 


{ 10 -6 ' 5 MgyU 1 t < 4 Myr 

10 ~ 5 ' 4 MgyU 1 ( 4 xl oe yr ) 4Myr < t < 40Myr. 

0 t> 40Myr. 

(B3) 

rhoB = fg.M/fhi,, where fg = 2.6 x 10 -3 is the fraction 
of the stellar mass with M* > 8 Mg and m* = 0.35 Mg is 
the mean stellar mass for our assumed Salpeter IMF, p ~ 
M*- 2 - 35 . 

For the mass loss rate from post-main sequence winds, 
we take 


m to 
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(B4) 


and 0 for earlier times. By truncating tHto at 40 Myr, 
we are ignoring the non-steady mass injection by type II 
supernovae. 

The quantity of mass lost in post-main sequence winds 
AM(f) is estimated from the expression given by Ciotti & 
Ostriker (2007) (their eq. [10]), 


AM 


0.945Mro - 0.503 Mro < 9M 0 
Mto — 1.4Mq Mro 7 9M 0 , 


(B5) 


where Mto is the turn-off mass, which at time t < th is 
calculated as 


log(Mro)/M 0 = 0.24 + 0.068a: 2 - 0.34a: + 4.76e“ 4 ' 58x , 

(B 6 ) 

where x = log(t/10 9 yr). This functional fit is designed to 
reproduce the results of Maeder & Meynet (1987) (their Ta¬ 
ble 9) for massive stars while asymptoting to the formula 
provided by Ciotti & Ostriker (2007) (their eq. [9]) for in¬ 
termediate and late times (t > 10 8 years). 

The first term in equation (B2) corresponds energy in¬ 
jection from massive stars, while the second term accounts 
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for energy injection from stars on the lower main sequence. 
Both terms have a thermalization efficiency, /, which w e 
talk to be 1. Note we do not account for energy from Type 
II supernovae. 

The MS wind mass loss rate rh(A/*,t) is calculated 
based on the generalization of Reimer’s law 

rh = 4x 10 " 13 ^ M ® yr " 1 ’ (B7) 

where R t , A*, T e g and <?* are the stellar radius, luminosity, 
effective temperature, and surface gravity, respectively, the 
latter normalized to its solar value gr 0 . The stellar radius and 
luminosity are estimated as (Kippenhahn & Weigert (1990); 
Figs. 22.2) 22.3) 


A* 


A 0 (A7*/M 0 ) 3 ' 2 Al* > M 0 
A 0 (A/*/M 0 ) 2 ' 5 Af* < M 0 


(B8) 


R* 


i? 0 (Af*/M 0 )°' 57 Af* > M 0 
R 0 (AA*/ M 0 )°' 8 Af* < M 0 


(B9) 


The wind velocity of main sequence winds is assumed to 
equal v w (M*,t) = u TO , 0 (Al*/M 0 ) 1 / 2 (f?*/Jf 0 ) -1 / 2 , i.e. scal¬ 
ing as the stellar escape velocity and normalized to the ve¬ 
locity of the solar wind v w ,@ = 430 km s _ ; this produces 
an effective wind heating velocity for main sequence winds 
alone of ~ 100 km s~ x for r* ~ th, close to the value found 
by Naiman et al. (2013) for globular clusters based on a 
more sophisticated population synthesis treatment. 

Winds from lower main sequence stars only dominate 
the energy injection a late times 10 Gyr after an impulsive 
burst of star formation. Even with 100% thermalization ef¬ 
ficiency these winds could not thermally stabilize the CNM. 

The rate of energy injection due to winds from massive 
stars, eoB(t) = fs£/fh*, where /g = 2.6xl0 ~ 3 is the fraction 
of the stellar mass with Af* > 8Af 0 for our assumed Salpeter 
IMF. Here £(t) is the energy injection rate per massive star, 
which we estimate as 




Figure Bl. Effective heating rate, u* , and mass loss parameter, 
?/ (eq. [8]), resulting from stellar winds from a stellar population 
of age r*. 


£(t) 


1.3 x 10 36 ergs _1 



t<4x 10 6 yr 

4 x 10 6 yr <t< 10 8 yrB2 Over Star Formation History 

t > 10 8 yr Generalizing to an arbitrary SFH S(t), the total rate of mass 

(B10) and energy input can be written as 


based on the results of Voss et al. (2009) (their Fig. 7, top 
panel), who use a population synthesis code to simulate the 
mass and energy injection into the ISM from an OB asso¬ 
ciation. Although equation is valid only for t < 10 Myr, in 
practice the precise functional form of eoB(f) is generally 
unimportant for our purposes, so we adopt equation as be¬ 
ing valid for all times. 

The effective wind velocity in the limit of an impulsive 
star formation may then be written as 

v w (t) = 2 d(t)/rh(t) (Bll) 

while 

r/ = fh(t)th (B12) 

Figure Bl shows the values of v w (t) and r](t) as a function 
of stellar age, r*. 


M(t) = f S(t\)m(t — fi)dti (B13) 

Jo 

E(t) — f S(ti)e(t — ti)dfi, (B14) 

Jo 

resulting in a wind heating parameter of 

vl(t) = 2 E{t)/M(t). (B15) 


The mass return parameter will be 


M(t) 

7) = —L - — - th 

Jo S(ti)dt i 


(B16) 


We estimate the stellar wind heating provided by the 
average star formation history of galaxies of a given M. us¬ 
ing the results of Moster et al. (2013, eqs. 17-20). Note that 
the star formation histories in Moster et al. (2013) are in 
terms of halo mass. For a given Af. we assign the halo mass 
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Figure B2. Effective wind velocities for nonstandard star forma¬ 
tion histories. The black curves shows, for reference, n w calculated 
using the halo-averaged S(t). and gray curves show wind heating 
resulting from perturbed star formation histories given by equa¬ 
tion (B24). In the top panel the star formation rate declines for 
a time, <5t* = 10 7 years to fractions t = 0.001 (dashed), i = 0.1 
(dotted), and l = 0.3 (dot-dashed) of the halo averaged value. 
The bottom panel shows results for <5t* = 10 s years and t=0.001 
(dashed) and 0.1 (dotted). 


whose star formation history would produce a bulge con¬ 
sistent with the M . — Mb u i ge relationship from McConnell 
& Ma (2013). A slight complication occurs for the largest 
mass halos, where much of the 2 = 0 stellar mass has been 
acquired through accretion of satellite halos rather than in 
situ star formation. To accommodate this, we incorporate 
analytic fits for mass accretion histories, taken from Moster 
et al. (2013, their eqs. 21-23), assuming that the age dis¬ 
tribution of the accreted stars is equal to the age distribu¬ 
tion of those formed in situ. This assumption may be con¬ 
servative if the primary galaxy’s accretion history is dom¬ 
inated by minor mergers with younger stellar populations. 
On the other hand, the dynamical friction inspiral time for 
small satellite galaxies is quite long, generally much greater 
than the ~ 10 7 yr for which young stars can dominate the 
heating budget. The mass of stars accreted for halo masses, 
Afhaio < 3 x 10 12 Mq, and redshifts, 2 > 4, is small and is 
neglected. 


To find the total mass ( M'(t )) and energy ( Eft )) injec¬ 
tion rates, including the contribution of accreted stars, we 
add a convolution of the specific mass and injection rates 
with the accretion history A(t) to the mass and energy in- 


jection rates from star formed in-situ. Thus, 


Mft ) 

= M(t ) + 

[ A(tacc)- 

Af(tacc) 

~ /. \ UL acc 

(B17) 

Jo 

f 0 S(ti)dti 

Eft) 

= m + j 

f A(t acc)- 

F'(iacc) j. 

UL acc 

(B18) 

'o J 

S(t')dti 


The corresponding wind heating parameter v' w and 
mass return parameter, r/, will be given by 


/ 

n 


M\t) 

fo S(ti)dti +fg A(ti)dti h 
_ 2 Eft) 

Mft) 


(B19) 

(B20) 


Figure B2 shows how the wind heating varies as star 
formation histories deviate from their halo-averaged values. 
In particular, we show 


M(t) = 

f S(ti)m(t — ti)dfi 

Jo 

(B21) 

m = 

[ - ti)dti, 

Jo 

(B22) 

V’(t) = 

2 £(t)/M{t) 

(B23) 


In these equations, 

S(t) = S(t) x ^ —(1 — i) arctan(f/5t*) + (B24) 

This function convolves the recent (z ~ 0) halo-averaged 
star formation history with local variation to give a more 
pessimistic estimate for the value of u w . In particular, re¬ 
placing S(t) with S(t) reduces the recent star formation rate 
to a fraction e of its halo-averaged value, and does so for a 
characteristic time 5t* into the past. As we can see in Fig. 
B2, this dramatically lowers the effective wind speed when 
both Sti, > 10 7 yr and t < 0.1, but otherwise has too modest 
of an effect to change the thermal stability properties of the 
flow (although the location of r s and the value of M may 
change significantly). 
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